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SUMMARY 


Mathematical models and an associated computer program have been devel- 
oped for heat pipe startup from the frozen state. The models have been checked 
against previously published analytical and experimental data. Agreement is rela- 
tively good for most situations examined. 

When a liquid metal heat pipe is started by introducing heat to one end while 
cooling the other, internal working fluid dynamics may greatly affect temperature 
distributions and fluid properties within the pipe as well as the overall conductance 
of the pipe. For example, if the working fluid is initially frozen, during startup 
melting will occur in the capillary structure and the vapor will experience free 
molecular, choked, and continuum flow at various times. These changing internal 
conditions generally make the heat pipe relatively slow to transport energy from 
heated to cooled ends and very large radial and axial temperature gradients may 
develop. 

The present work uses finite element formulations of the governing equations 
written for each heat pipe region for each operating condition experienced during 
startup from a frozen state. In the shell, energy transport is by conduction only. In 
the capillary structure, conduction and heat of fusion are considered. In the vapor 
region different sets of governing equations axe utilized for regions undergoing free 
molecular, choked and normal continuum flow. The various models were checked 
against analytical and experimental data available in the literature for three specific 
types of operation. For example the models used to predict melting in the capillary 
structure were checked against analytical results previously published for melting 
in a corner region. 

Computation using the methods developed in the present work were made for 



a space shuttle reentry mission where a heat pipe cooled leading edge was used on 
the wing. This wing had a sodium heat pipe built into the wing near the leading 
edge. Charles J. Camarda of NASA Langley Research Center made experimental 
measurements of startup behavior for such a heat pipe. Results computed in this 
work compared well with Camarda’s data. 
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CHAPTER I 
INTRODUCTION 

1.1 Statement of the problem 

Since the heat pipe concept was first introduced independently by Gaugler[l] 
in 1942 and Grover[2] in 1963, theoretical and experimental studies have been un- 
derway to understand and develop heat pipe technology. Most theoretical stud- 
ies concern certain portions of the heat pipe, such as the evaporator, condenser, 
capillary structure, and vapor flow region. The overall performance of the entire 
heat pipe, including the thermal behavior along the heat pipe wall and capillary 
structure, vapor flow dynamics, and the various types of boundary conditions on 
the evaporator and condenser surfaces have received less attention. However, the 
steady state characteristics of heat pipe performance at low temperatures and un- 
der normal operating conditions are relatively well understood, and heat pipes have 
been successfully applied in various fields. 

Little research has been done on the transient case. Transient behavior of heat 
pipes have been experimentally and numerically studied for low temperatures and 
working fluids with high vapor density by Chang and Col well [3 ,4, 5]. 

Recently, use of the heat pipe has been considered as a means of reducing the 
peak temperature and alleviating the thermal gradients at the leading edges of reen- 
try vehicles and hypersonic aircraft, and in nuclear reactors. In these applications, 
the rate of heat transfer may be large, and the range of operating temperatures 
broad, from ambient to high temperature, so that liquid metal, which is in the solid 
state at ambient temperature, may be used as the working fluid. Under these condi- 
tions, the working fluid in the capillary structure may be in the solid or liquid state, 
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or may be freezing or thawing, with some liquid and some solid present. The vapor 
flow may be free molecular, continuum, choked, or some combination of these. No 
complete research has been reported on this subject. 

The primary objective of this research is to investigate analytically the per- 
formance of an entire heat pipe with metallic working fluid during startup from a 
frozen state. To accomplish this goal, a mathematical model has been developed, 
and a numerical solution technique tested to predict the transient temperature dis- 
tributions along the heat pipe, and the optimal heat transfer rate. 

1.2 Description of heat pipe operation 

Many scientists and engineers have observed the phenomenon of surface tension 
in nature and tried to understand, formulate, and apply it for improving human 
life. Among many natural phenomena, the action of surface tension can raise liquid 
against a gravity force within a small vertical tube or gauze with a portion immersed 
in liquid. This capillary action can transport liquid through suitable materials with- 
out using external power. When phase change takes place from one state to another, 
the change of enthalpy is rapid, and the difference between enthalpies of two states 
is large. Therefore a large amount of energy may be absorbed or released depending 
on the direction of phase change, and without a large temperature gradient. These 
two phenomena are utilized in a heat pipe. 

A heat pipe consists mainly of a shell as the container, a capillary structure or 
wick to transport liquid by using surface tension, and a vapor space to provide vapor 
passage as shown in Figure 1.1. Heat pipe shells have been made of stainless steel, 
copper, nickel alloys, hastelloy, et cetera. Wire screen, fiber glass, porous metal, and 
w r oven cloth have been used as capillary structures. Narrow' grooves cut lengthwise 
in the interior pipe wall have also served as a capillary structure. The capillary 
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structure is saturated by the working fluid in the liquid state, and the vapor space 
is occupied by the working fluid in the vapor state. Heating and/or cooling devices 
axe applied to the outer surface of the heat pipe shell. In the longitudinal direction, 
the heat pipe consists of an evaporator and a condenser. 

Heat added to the evaporator is transferred to the working fluid by conduction 
and causes vaporization of the working fluid at the surface of the capillary structure. 
Vaporization causes the local vapor pressure in the evaporator to increase and vapor 
to flow towards the condenser thereby transporting the latent heat of vaporization. 
Since energy is extracted at the condenser, the vapor transported through the vapor 
space is condensed at the surface of the capillary structure, releasing the latent 
heat. The radius of curvature of the meniscus in the capillary structure of the 
evaporator is decreased and that in the condenser is increased. This difference in 
radii between the two sections creates the pumping force that transports the liquid 
from the condenser to the evaporator through the capillary structure. This process 
continues so long as no extreme heat fluxes axe encountered. 

Hence, in a heat pipe energy is transported by utilizing phase change of the 
working substance instead of a large temperature gradient and without external 
power. Also, the amount of energy transferred through a small cross-section is 
much larger than that by conduction or convection. Heat pipes may be operated 
over a broad range of temperatures by choosing an appropriate working fluid, as 
shown in Table 1.1 [6]. 

However, this useful device has some operating limitations such as the sonic 
limit, the capillary Emit, the entrainment limit, and the boiling Emit. When any 
of these Emitations is encountered, the capillary structure may dry out leading to 
failure of the heat pipe. In addition to these Emitions, when Equid metal is used 
as the working fluid, startup difficulty may take place due to possible solid state of 
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Table 1.1 Heat pipe working fluids from Ref.[6] 


Medium 

Melting 

point 

(°C) 

Boiling point 
at atmos. press. 
(°C) 

Useful range 
(°C) 

Helium 

- 272 

- 269 

- 271 - 

- 

269 

Nitrogen 

- 210 

- 196 

- 203 - 

- 

160 

Ammonia 

- 78 

- 33 

- 60 


100 

Freon 11 

- Ill 

24 

- 40 

- 

120 

Pentane 

- 130 

28 

- 20 


120 

Freon 113 

- 35 

48 

- 10 

- 

100 

Acetone 

- 95 

57 

0 

- 

120 

Methanol 

- 98 

64 

10 

- 

130 

Ethanol 

- 112 

78 

0 

- 

130 

Heptane 

- 90 

98 

0 


150 

Water 

0 

100 

30 

- 

200 

Toluene 

- 95 

110 

50 

- 

200 

Mercury 

- 39 

361 

250 

- 

650 

Cesium 

29 

670 

450 

- 

900 

Potassium 

62 

774 

500 

- 

1000 

Sodium 

98 

892 

600 

- 

1200 

Lithium 

179 

1340 

1000 

- 

1800 
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the working fluid and extremely low vapor density. 

1.3 Literature review 

The heat pipe is a highly effective device for transporting heat between a source 
and a sink. Since Gaugler[l] received a patent on the heat pipe concept applied to 
a refrigeration system and Grover[2] referred to a “heat pipe” in a patent filed for 
the United States Atomic Energy Commission in 1963, scientists and engineers have 
been developing heat pipe technology. The first paper which described the basic 
principle of operation of a heat pipe was published by scientists at the Los Alamos 
National Laboratory[7] in 1964. They built two heat pipes with water and sodium 
as the working fluids for an initial qualitative experiment. Work at Los Alamos[8,9] 
continued actively, emphasizing space applications for the transfer of very high heat 
fluxes between two components and for the elimination of temperature gradients 
over relatively large areas. For high temperature applications, lithium and silver 
were tested as the working fluid at 1300° C and 2000° C, respectively. For the first 
actual flight test, a stainless steel heat pipe with distilled water as the working fluid 
operated successfully. At this stage, research on heat pipes was also conducted in 
England and Italy [6]. 

Cotter[10] developed the general basic theory for making certain quantitative 
calculations of heat pipe behavior. This analysis was confined to right circular 
cylinders of large length-to-diameter ratio and emphasized the vapor flow. Uniform 
injection or suction were assumed for a steady state condition. The axial transport 
of energy was modeled with the vapor flow carrying the latent heat of vaporization 
while neglecting axial conduction and radiation in the vapor space. After this 
pioneering effort, several books[6,ll,12] were published which describe the basic 
theory of conventional heat pipe operation at steady state and low temperature. 
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Busse[13] studied the pressure drop in laminar vapor flow in a long, cylindri- 
cal heat pipe. The vapor density was assumed to be constant. The Navier-Stokes 
equations were simplified by a boundary layer approximation and solved by approx- 
imating the axial velocity component as a polynomial of the fourth power of the 
radius, with a correction function which varied only in the flow direction. 

Levy[14,15] used a one-dimensional approach in solving the vapor flow problem, 
taking into account compressibility, shear stress at the liquid-vapor interface and 
the vapor dissociation-recombination reaction. The analytical results indicated that 
the shear stress is the most important factor, which reduces the maximum rate of 
heat transfer from that based on the sonic limit. 

Brovalsky et al.[16] described the vapor flow for alkali-metals by using averaged 
equations of motion over the cross-section. Compressibility and friction at the 
liquid-vapor interface were considered. Momentum and energy factors, and the 
friction factor were evaluated based on theoretical data available for incompressible 
vapor flow in a channel with porous walls. The comparison of numerical results 
with available experimental data indicated a maximum discrepancy of 10 %. The 
temperature drop along the heat pipe was also observed. 

Bankston, and Smith[17] studied the fluid dynamics of the vapor flow at three 
different Reynolds numbers; 0.01, 4, and 1000. The Navier-Stokes equations for 
steady, incompressible, laminar vapor flow in a cylindrical heat pipe were solved 
by a finite difference method in which the dependent variables were transformed 
to the stream function and the vorticity. Inflow and outflow boundary conditions 
were described at the wall as blowing and suction through a porous wall pipe, but 
no thermodynamic change of phase was actually employed in the analysis. Their 
results show that vapor flow’ in the condenser is more complex than that in the 
evaporator, and that the vapor pressure varies not only in the axial, but also in the 
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radial direction for large Reynolds numbers. 

A numerical analysis of steady two dimensional heat and mass transfer in the 
vapor-gas region of a gas loaded heat pipe was made by Tien and Rohani[18]. In this 
study, the radial component of velocity at the wall was determined by writing an 
energy balance equation at the liquid-vapor interface. Numerical results show that 
considerable pressure variations in the axial direction exist for large heat fluxes, 
due to friction at the liquid-vapor interface. Thus, a temperature drop in the 
axial direction occurs and vapor pressure variations play a significant role in overall 
performance of the heat pipe. 

Vapor flow dynamics in a flat plate heat pipe with asymmetric boundary con- 
ditions was studied by Ooijen and Hoogendoorn[19]. The numerical study was con- 
fined to two-dimensional, steady state, laminar and incompressible flow. From com- 
putational results, velocity profiles were plotted in the evaporator and condenser for 
wall Reynolds numbers of 2, 10, and 50, and were compared with parabolic Poiseuille 
profiles for three locations. Flow reversal was observed, and similarity did not ex- 
ist for high wall Reynolds numbers in the condenser section. For small Reynolds 
numbers, the pressure drop is similar to the Poiseuille flow model. However, for 
high wall Reynolds number, the large difference in velocity profiles in the condenser 
section causes a higher pressure drop than that resulting from the Poiseuille model. 
Good agreement was observed between experimental and computational results for 
nitrogen gas. 

Ismail and Murcia[20] studied combined liquid and vapor flow in a tube with a 
porous wick. Governing equations for the flow of viscous incompressible fluid were 
solved using the separation of variables with known evaporation or condensation 
rates. For the case of small Reynolds number, analytical results were obtained. 

Demichele[21] investigated the two-dimensional, steady state and compressible 
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flow problem by using an integral transformation of the general compressible flow 
equation introducing stream tubes which can be thought of as a set of concentric 
nozzles. For each stream tube equations were derived with a different initial con- 
dition. Numerical solutions predicted velocity, pressure, and temperature profiles. 
Effects of viscous terms on pressure recovery were deemed to be important. 

Compressible, transient and axisymmetric Navier-Stokes equations were nu- 
merically solved to derive a friction coefficient expression to be used in one-dimens- 
ional heat pipe vapor models by Bowman[22]. The equation for the friction coef- 
ficient was expressed in terms of local axial Reynolds number, Mach number, pipe 
aspect ratio and radial Reynolds number. The one- dimensional model with the 
friction coefficient showed good agreement with experimental results. 

Overall thermal performance of a heat pipe at steady state was studied the- 
oretically and experimentally, by Sun and Tien[23,24]. In the analysis, a simple 
conduction model was developed for a single-component heat pipe in one dimension 
with uniform saturated vapor temperature and uniform mass injection or suction. 
Axial wall temperature distributions were predicted. Theoretical predictions were 
compared with measured results for gas-loaded heat pipes and good agreement was 
reached. 

As the digital computer was developed numerical techniques were used to solve 
more complicated models. The simple conduction model developed by Sun and Tien 
was extended by Kuramae[25] to transient heat pipes with time varying thermal 
loads. In this model, the temperature was assumed to vary only in the axial direction 
for the wall but in both the axial and radial directions for the wick structure. It was 
assumed that the vapor temperature was dependent on time but uniform in space. 
A numerical method was used to solve the governing equations and the calculated 
axial temperatures were compared with typical experimental results. 
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Experimental and numerical studies were conducted for transient operating 
characteristics of low temperature heat pipes by Chang and Colwell[3,4,5]. The 
computational model assumed that two-dimensional conduction in the heat pipe 
shell and wick are the dominant heat transfer modes. A lumped mass model was 
used for a combination of vapor space and central composite slab wick. Thermal 
resistance at the liquid-vapor interface and along the vapor space was neglected. 
To provide boundary conditions at the interface between the heat pipe wall and the 
vapor space, the vapor region was modeled assuming that the vapor temperature was 
a function of time only. A finite-difference method was used to predict performance. 

Cotter[26] described three basic transient modes for heat pipe startup. A 
frontal startup mode was observed when the vapor density is so low that the molec- 
ular mean free path exceeds the diameter of the vapor passage. In this mode of 
startup, the vapor in the hot zone is in continuum flow and that in the cold zone 
is in free molecule flow. A large temperature gradient is developed and decreases 
with time. Eventually, an isothermal steady state could be reached. 

Ivanovskii et al.[ll] carried out experimental studies of the temperature distri- 
bution along the length of a sodium heat pipe in which the working fluid was in the 
solid state at ambient temperature. The temperature distributions were measured 
with the aid of a movable microthermocouple placed directly in the vapor channel. 
Ivanovskii observed three simultaneous flow regimes in the condensation zone for 
intense heat removal in a pipe operating at the sonic heat transfer limit: continuum 
vapor flow at the start of the condensation zone and intermediate and free molecular 
regimes futher on. 

Neal[27] investigated the successful startup of a heat pipe with water as the 
working fluid. The water in the condenser was initially frozen, but that in the 
evaporator remained in the liquid phase. Experimental results showed that a large 
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temperature gradient developed along the heat pipe length with increasing time 
and even for a small heat flux the heat pipe failed to start up. 

Shlosinger[28] studied the startup behavior of low temperature heat pipes with 
the water initially frozen. With a heat input of 15 watts, the heat pipe working 
fluid melted without wick dryout and normal operation began after approximately 
one hour. Use of an auxilliary heat pipe which had a working fluid with a lower 
melting point greatly improved the startup of the primary heat pipe without local 
overheating and the transient period was reduced. 

Deverall et al.[29] made a series of tests to study the startup problem with water 
and metallic heat pipes. They described the transient behavior of heat pipe startup 
based on their experiments. With the working fluid in the solid state, startup was 
possible, but was highly dependent on the heat rejection rate at the condenser. For 
successful startup, the heat rejection rate at the condenser had to be low enough 
to enable the heat to melt the working fluid in the condenser, and allow liquid to 
return through the wick structure before the evaporator was depleted of fluid. Heat 
rejection by radiation is a self-compensating system and automatically controls the 
heat rejection rate. Therefore, startup difficulties were not normally encountered 
under these conditions. Another method suggested to aid the startup of a heat pipe 
was the addition of a small amount of inert noncondensable gas which has a result 
similar to that of radiation. 

Camarda[30] investigated the performance of a heat pipe cooled leading edge, 
experimentally and analytically. In the analysis, it was assumed that the temper- 
ature was uniform throughout the continuum flow region and the melting process 
of the working fluid was neglected. Temperatures were calculated by a lumped sys- 
tem method which used a volumetric heat capacity per unit length of heat pipe. 
Rates of continuum vapor region growth, which were predicted using simple energy 
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balances, were compared with experimental results. 

Most experimental studies[27,28,29] attempted simply to see if it was possible 
to obtain successful startup from the frozen state. It has been noted that startup 
difficulties are normally encountered when the working fluid is initially in the solid 
phase. No analytical studies which include the effects of phase change of the working 
fluid and vapor flow dynamics have been reported, and further more, comprehensive 
and qualitative research on startup from the frozen state has not been carried out 
experimentally and numerically. 
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CHAPTER II 

MODELING OF STARTUP 

2.1 Introduction 

The use of heat pipes is being considered as a means of reducing the peak 
temperature and large thermal gradients at leading edges of reentry vehicles and 
hypersonic aircraft, and in nuclear reactors[31,32]. In the basic cooling concept, 
the heat pipe, which is a highly effective heat transport device, covers the leading 
edge, a portion of the lower wing surfaces, and a portion of the upper wing surface. 
Aerodynamic heat is mainly absorbed at the leading edge, and transported through 
the heat pipe to the upper and lower wing surfaces, where it is rejected by thermal 
radiation and convection. Once fully operational, the near isothermal heat pipe 
virtually eliminates temperature gradients and reduces peak temperatures. 

A previous feasibility study[31] of heat pipes for this application recommends 
a rectangular cross-section for the heat pipe based on weight per unit surface area 
of heat pipe cooling structure and fabrication considerations. A schematic diagram 
of the physical model based on results presented in reference[31] is shown in Figure 
2 . 1 . 

Previous experimental observations[26-32] suggest the following sequence of 
events during heat pipe startup from the frozen sate. Initially, the working fluid is 
in the solid state and the vapor density is extremely low, so that free molecular flow 
conditions prevail throughout the vapor space. The input flux over the evaporator 
starts to melt the frozen substance in this region, while the heat transport from the 
hot zone to the adjacent zone proceeds quite slowly via axial conduction through 
working fluid and capillary structure, while heat transfer in the vapor is almost 
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negligible. 

When energy is continuously added to the evaporator, the frozen working fluid 
in the evaporator is melted, so that evaporation can take place at the liquid-vapor 
interface and vapor flows into the condenser section due to the large pressure gra- 
dient. Vapor therefore freezes on the inner surface of the frozen working fluid in 
the cold zone and the vapor-solid interface temperature increases until the melting 
temperature is reached. During this stage, energy is mainly transferred as latent 
heat owing to vaporization in the heated zone, and condensation and freezing in the 
cooled zone. The vapor flow may be choked at the exit of the evaporator because 
of very low pressure in the cold zone. 

This process continues until the frozen working fluid is completely melted and 
the continuum flow regime reaches the end of the heat pipe, at which time liquid 
returned to the evaporator is sufficient for normal transient operation. Eventually 
the heat pipe may reach a steady state condition. As suggested, during the startup 
of the heat pipe from a frozen state, the behavior of vapor flow may be divided into 
three distinct phases for convience of analysis. 

Phase Is Vapor flow in the heat pipe is in free molecular condition through the 
vapor space. 

Phase II: In the vapor space, a region of continuum flow is established in the 

heated zone and a continuum flow front moves toward the heat pipe 
cooled end. Vapor flow may be choked at the end of the evaporator. 
Phase III: Continuum flow exists over the entire heat pipe length in the vapor 
region and the sonic limit is not encountered. 



16 


2.2 Mathematical model development 

On the basis of experimental observations, basic governing equations are writ- 
ten to determine the startup, transient, and steady state performance of a heat pipe 
which has initially frozen alkali-metal as the working fluid. These equations can be 
coupled by several types of boundary conditions on the heat pipe surface, such 
as specified temperature, heat flux, convection and radiation boundary conditions. 
The boundary condition at the liquid-vapor interface depends on the three phases 
of vapor flow dynamics mentioned in section 2.1. 

2.2.1 Transition temperature 

Continuum flow in the vapor space is considered to be established when the 
mean free path, A, is substantially less than the minimum dimension,!?, of the vapor 
flow passage, e.g., 


K n = -< 0.01 


( 2 . 1 ) 


From kinetic theory of gases[33], the dynamic viscosity and the mean molecular 
velocity can be expressed as 


fjL = 0.5pAF 


( 2 . 2 ) 


V = 



(2.3) 


Eliminating the mean free path from Equations(2.1) and (2.2) yields the tem- 


perature of the vapor space corresponding to the given mean free path, 
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T* > 


7T M 
2 x 10-4^ 



(2.4) 


where p is the density of the vapor, p is the dynamic viscosity of the vapor, R u is 
the universal gas constant, and M is the molecule weight. 

Iterations are required to obtain a value of T* due to the temperature depen- 
dence of properties. Figure 2.2 illustrates the transition temperature, T*, from free 
molecular to continuum flow as a function of minimum vapor passage for sodium. 
When the temperature of the vapor space is greater than that calculated by Equa- 
tion(2.4), continuum flow is assumed to be established in the vapor space. 

2.2.2 Heat pine shell and capillary structure with working fluid 

Unlike the case of a conventional cylindrical heat pipe, aerodynamic heating 
for the heat pipe shown in Figure 2.1 will cause the highest temperature to occur at 
the outer shell of the heat pipe. Thus, to simplify the analysis, it is assumed that 
the inner shell is sufficiently thin to neglect its thermal resistance and capacitance. 
Furthermore, the following additional assumptions axe made. The thicknesses of 
the heat pipe and the rib are assumed to be much smaller than the width and 
therefore temperature gradients would be developed primarily in the chordwise 
direction during heat pipe startup. Therefore, the rib may be neglected in the 
analysis. Also, the heat pipe is assumed to be symmetric about the stagnation line 
so that the upper half section of the heat pipe only is considered. Hence, a two- 
dimensional model is considered for mathematical formulation. A simplified cross 
section is shown in Figure 2.3. 

The unsteady, two-dimensional conduction equation is applied to the heat pipe 
shell and the capillary structure to evaluate the temperature distributions under 
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Outer Wall 


Liquid and Wick 


Vapor Space 


Adiabatic boundary condition 


Figure 2.3. Simplified cross-section A- A of Figure 2.1. 
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the following assumptions: 

1. The heat transferred through the wick and working fluid is by conduction 
only since liquid flow velocity is very low, liquid layer is thin and thermal 
conductivity is very high; 

2. The presence of the wick structure does not affect the melting of the work- 
ing substance, and melting occurs over a small, finite temperature range, 
AT, around the melting temperature; 

3. Boundary conditions at the liquid-vapor interface are subject to the phases 
of vapor flow dynamics noted in section 2.1. 

4. Radiation heat transfer at the liquid-vapor interface is neglected based on 
the small thermal emissivity of liquid sodium[34]. 

The governing equations and boundary conditions for the heat pipe shell and 
capillary structure with working fluid are expressed in one form: 






i = 1,2 


(2.5) 


T{ = T 0 for t < 0 


( 2 . 6 ) 


K : 


m 

dn 


Q(X,t ) at Y = fi(X) for 0 < i/> < V> e 


(2.7) 


-K i 


dTj 

dn 


k C r(T-i - T cr ) + cre(Ti 4 - T r 4 ) at Y = f\(X) for < A 

( 2 . 8 ) 


Ti 


T 2 


and 


Ki 


dT\ 

dn 


= K 2 


dT 2 

dn 


at Y = f 2 (X) for 0 < -0 < i' c (2.9) 
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dTi 

ds 


= o 


i = 1,2 


at ip = 0 and ip = ip c 


( 2 . 10 ) 


where the subscripts 1 and 2 denote the heat pipe shell and the capillary structure, 
n represents the unit outward normal direction, and Sd represents the chordwise 
direction. 

In X - Y coordinates 


ar ar 

K dn K ax' 


+ IC ?L 

T A ar " 


( 2 . 11 ) 


where l x and l y are the direction cosines between the surface outward normal, n, 
and the X and Y axes, respectively. 


2.2.3 Analysis of the vapor flow 


- Analysis of the vapor flow is necessary to provide the boundary condition at 

the liquid-vapor interface when continuum flow is established in part of the vapor 
space. However, the behavior of the vapor flow as mentioned in section 2.1 is very 
complicated due to the extremely small vapor density of the metallic working fluid 
at low vapor pressure, and the large pressure gradient in the chordwise direction. 
Limits on vapor flow are encountered, and considerable thermal resistance exists at 
the liquid-vapor interface due to phase change of the working fluid. The case of a 
flat plate heat pipe with asymmetrical boundary conditions shown in Figure 2.3 is 
considered for the vapor flow passage. 

2. 2. 3.1 Evaporation and condensation 


From kinetic theory, a significant thermal resistance exists at the liquid-vapor 
interface for liquid metal[35] and this resistance increases with decreasing vapor 
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pressure. This means that the difference between the vapor temperature and the 
interface temperature is not negligible. At the interface, interphase mass transfer 
has been stated from the viewpoint of kinetic theory as a net molecular flux which is 
the difference between the rate of arrival of molecules from the vapor space towards 
the interface and vice-versa. While condensation is proceeding, the arrival rate of 
molecules exceeds the departure rate. During evaporation, the reverse is true, and 
during the equilibrium state, the net molecular flux is zero. Hence, evaporation and 
condensation are modeled by using an expression for the net rate derived from the 
kinetic theory of gases[36]: 


rh 0 


( 2*0 / M 

j 

\2 — a ) V 27r.Ru 

>/Tf y/T g \ 


( 2 . 12 ) 


where a is the condensation or evaporation coefficient which is assumed to be unity, 
e is the porosity of the wick, rh 0 is the rate of condensation or evaporation per unit 
area, M is the molecular weight, R u is the universal gas constant, Pf and 7/ sure 
the pressure and temperature, respectively, at the interface, and P g and T g are the 
pressure and temperature of the vapor, respectively. 


2. 2. 3.2 Limitations of vapor flow 

After continuum flow is established in the vapor space, because of the low 
density of vapor at low pressure and the large pressure gradient, the vapor flow is 
choked at the end of the evaporator, even for a low heat flux, just as it is at the 
throat of a convergent nozzle for large pressure gradient. This sonic limit is the 
first among several limitations encountered, and the performance of the heat pipe 
is restricted until the vapor temperature increases accordingly until the velocity of 
the vapor leaving the evaporator is less than the sonic velocity. 
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The expression for the sonic limit[15,29] in terms of flow conditions at the 
beginning of the evaporator, from an energy balance on a control volume enclosing 
the entire evaporator, is: 



(2.13) 


where A c is the cross-section area of the vapor space, hf g is the enthalpy of vapor- 
ization, 7 is the ratio of specific heats, M is the weight of a molecule, R u is the 
universal gas constant, T 0 is the vapor temperature at the beginning of the evapo- 
rator and p 0 is taken as the saturation density corresponding to T 0 . Experimental 
verification of Equation(2.13) has been made by Kemme[37] for sodium, potassium 
and cesium heat pipes. 

Figure 2.4 shows the sonic limit, the entrainment limit[12} and the axial Rey- 
nolds number for a heat pipe which has a cross-sectional area of 0.55 cm 2 , hydraulic 
radius of 0.32 cm, 100 mesh screen wick, and sodium as the working fluid. The type 
of limitation restricting performance of a heat pipe is determined by which limitation 
has the lowest value at the temperature under consideration. Thus, the first limit 
encountered is the sonic limit, as shown in Figure 2.4. However, when the actual 
chordwise heat transfer required is below these limits, no limits are encountered. 


2. 2. 3. 3 Modeline of vapor flow durine phase II 


Even though some of the working substance is in the liquid state, the transition 
temperature from free molecular flow to continuum flow is much greater than the 
melting temperature, so that the vapor flow is assumed to be free molecular during 
phase I, and an adiabatic boundary condition is applicable at the interface. How- 
ever, during phase II, a region of continuum flow is assumed to be established in the 
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adjacent vapor space when the interface temperature is greater than the transition 
temperature, while in the cold zone the vapor is still in free molecular flow. The 
continuum flow region grows until it reaches the end of the heat pipe. The temper- 
ature in the region of free molecular flow remains unchanged except in the vicinity 
of the continuum flow region. Therefore, an imaginary plane, which is adiabatic 
and normal to the liquid-vapor interface, is assumed to divide the two vapor flow 
regions at the point of the transition temperature, and the dividing plane moves 
toward the cooled end of the heat pipe as the location of the transition temperature 
at the interface moves. 

In the continuum flow region, energy is mainly transported as latent heat of the 
working fluid, and variations of temperature and pressure in the chordwise direction 
are significant. Even though continuum flow exists in the vapor space, the vapor 
pressure is low and the pressure gradient in the chordwise direction is large, so the 
heat transfer is limited by the choked flow condition, and supersonic vapor flow and 
a shock front [6, 11, 23, 29] may occur in the condenser. 

In this research, the overall performance of a heat pipe is of more interest 
than details of one portion, and as noted, the ultimate heat transfer rate in the 
axial direction is limited by the sonic limit. Thus, the vapor flow during phase II 
may be approximated, without losing generality and accuracy, by evaluating the 
sonic limit properly. In order to evaluate the sonic limit, the total heat input to 
the vapor space, which can be obtained by applying Equation(2.12) to the element 
adjacent to the continuum flow region, is equated to the sonic limited transport 
from Equation(2.13) as follows: 


Tfl e 

I> 


/ M 
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/ 2t tR u 
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V A R, t T, 


(2.14) 
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where A c is the cross-section area, A Li is the length of the element, m e is the 
number of elements at the interface in the evaporator, and W is the width of the 
vapor space. From Equation(2.14), the uniform vapor temperature can be computed 
by iteration. This vapor temperature may be lower than true vapor temperature 
at the beginning of the evaporator, but at this stage the variation of density with 
temperature and the density itself is very small, so that the variation of sonic limit 
with respect to temperature is also small, as shown in Figure 2.4. Thus, the heat 
flux at the liquid-vapor interface and the sonic limit can be obtained by solving 
Equations(2.12) and (2.13), respectively, with the vapor temperature obtained from 
Equation(2.14). This method may be used until the vapor flow state reaches phase 
III. 


2.2.3.4 Vanor flow during phase III 

In this stage, the entire working fluid is completely melted and continuum 
flow exists throughout the vapor space. However, the heat pipe has not reached 
the desired operating condition and the density of the vapor is still small, so that 
compressibility should be considered. The amount of energy stored in the vapor 
space is negligible. The Reynolds number in Figure 2.4 is the maximum value 
corresponding to a given temperature and the diameter of the vapor space. When 
the actual heat transfer required is below these limits, no limits are encountered. 
A model[30] test, in which a thermal history for space shuttle reentry heating was 
simulated, showed that the typical maximum heat pipe operating temperature is 
about 940 K. At this temperature, the maximum Reynolds number is 1200. A study 
of the effect of mass injection and suction, and/or chordwise pressure gradients on 
flow transition from laminar to turbulent flow in the vapor space of a heat pipe, 
was conducted by Bowman[22]. His results showed that mass injection and pressure 
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drops in the evaporator correspond to those for laminar flow for axial Reynolds 
number up to 10 6 . For the condenser, transition of flow is observed at a Reynolds 
number of 12,000. Even though the Reynolds number depends on the geometry of 
the heat pipe and the actual heat transfer rate, the results cited from the previous 
study show that the vapor flow may be assumed laminar. Also, it was observed 
that the vapor reaches the steady state quickly, while the thermal response of the 
heat pipe wall and wick progresses slowly. It was recommended that a steady state 
model can be employed for the vapor flow. 

Thus, quasi-steady, compressible, one-dimensional laminar flow in the vapor 
space is considered. For purposes of formulating the mass, momentum, and energy 
equations in one-dimensional form, the velocity is taken to be the average velocity, 
which is approximated by the velocity distributions based on the similarity solutions 
of semiporous channels. In addition, friction at the liquid-vapor interface, variations 
of vapor quality, and momentum and energy factors axe similarly calculated. The 
vapor is evaporated at the interface with mass injection rate per unit area,m 0 . It 
is assumed that this vapor flows inward with a normal component of velocity only, 
and joins the chordwise vapor flow. 

Mass, momentum and energy balances in the chordwise direction, with the 
assumptions noted, yield: 


D jg(p v ) = ™° 


(2.15) 


dP d . .. 
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where D is the height of the vapor space, h is the vapor enthalpy and h 0 is the 
vapor enthalpy at the interface. 

A friction factor F for the surface is written as 


_ 8 Tg 

F = 9 


pV 2 


( 2 . 18 ) 


Momentum and energy factors, Mf and Ef, respectively, are expressed as 


M f 


[ U 2 dy 
J 

(2.19) 
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fk 


U 3 dy 

(2.20) 


Normal velocity of the vapor at the interface is expressed in terms of the heat 
transfer rate and latent heat of vaporization as follows: 


where Q is the heat 
temperature and A 0 


V 0 = 


Q 

hfgPo-Ao 


( 2 . 21 ) 


input rate at the interface, p 0 is the density for the interface 
is the interface area. 


i 


29 


CHAPTER III 

NUMERICAL MODEL DEVELOPMENT FOR PHASE CHANGE 

3.1 Introduction 

Transient and nonlinear heat transfer problems having phase change have been 
encountered in many engineering fields, and pose inherent difficulties for analytical 
and numerical solutions, because the interface between the new and old phase is 
moving with time, and properties such as conductivity, specific heat, and density 
are discontinuous at the phase change region. 

The fundamental feature of these phase change' problems was given attention 
and was solved analytically by Lame and Clapeyron in 1831, Stefan in 1891, and 
Newmann in 1912. Since then, many scientists have introduced analytical methods 
of solution of phase change problems to a number of idealized problems involving 
semi-infinite or infinite regions, subject to simple boundary and initial conditions. 
A brief discussion of these methods is presented by Ozisik[38]. A large number of 
numerical solutions of phase change problems were made possible by the availability 
of high speed digital computers. The finite element method has been used for 
nonstructural problems since the procedure was first proposed by Zienkiewicz and 
Cheung[39]. The solution process is now well established for linear situations, but 
relatively little work has been reported for nonlinear problems. 

The Galerkin weighted residual method is used here to derive finite element 
formulations. Since the governing differential equations are highly nonlinear due 
to the temperature dependence of the thermophysical properties, a three time level 
difference scheme which was proposed by Dupont et al.[40] is utilized to allow a 
direct evaluation of the properties at an intermediate time level, thereby eliminating 
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iterations within each time step. However, this scheme is not self-starting, so the 
implicit method with Newton-Raphson iteration is used for the first few time steps. 


3.2 Mathematical model 


For the system, having am arbitrary control volume V and control surface area 
A, the general principle of conservation of energy implies that whatever energy 
enters the system must either leave the system across the boundary or cause a 
change in the energy within the system. With no work, mass transfer, nor energy 
generation sources, the net rate of change of the total energy within the system 
must equal to the net rate of energy entering the system due to heat transfer across 
the control surface area. This statement can be expressed in mathematical form as 
follows: 


^ [ phdV = f KVT • ndA (3.1) 

di Jv Ja 

where p is the density, h is the specific enthalpy, K is the thermal conductivity, VT 
is the gradient of temperature and n is the unit outward normal to the surface as 
shown in Figure 3.1. 

For a single phase region, the fixed control volume is independent of time, and 
the divergence theorem is used to convert the surface integral to a volume integral. 
Therefore, Equation(3.1) can be written 


Ui {ph)dV= L 


V • (KVT)dV 


(3.2) 


The volume may be chosen so small as to remove the integral, and the specific 
enthalpy can be replaced by the expression for the specific heat and temperature. 
Equation(3.2) is then reduced for an arbitrary small volume V as follows: 
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Figure 3.1. Control volume V'. 
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pCp ^± = V . (KVT) (3.3) 

This equation can be applied to the solid and liquid regions when the motion 
of the liquid due to the change of the density is neglected. However, the properties 
of the material are discontinuous at the interface as shown in Figure 3.2 , so that a 
single equation cannot describe the phase change phenomenon. Hence, a problem 
with the phase change of a substance is mathematically described as follows: 


C21 ^ = V • (K21 vr 21 ) for old phase (3.4) 

at 

flrp 

C 22 — = V • (K 22 VI 22 ) for new phase (3.5) 

ot 

with initial condition 

T n = To (3.6) 

Boundary conditions are, for i = 1,2, 

T 2i = T„ on A\ (3.7) 

-K 2i VT 2i = Q on A 2 (3.8) 

-K 2i VT 2i = h cr (T 2 i - T cr ) on 4* (3.9) 


-K 2 i'VT 2 i — /3 r (T 2 i — T r ) on A 4 


(3.10) 
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where = ae(T 2 i 2 + T r 2 )(T 2 i + T r ) and coupling conditions at the interface axe 


T 2 i = T 22 = T m (3.11) 


K dT 2 1 K dT 2 1 _ , T dSi , . 

where h cr is the heat transfer coefficient, Hi, is the latent heat of a phase change, 
S is the interface position, T m is the temperature of the phase change, t is the 
emissivity, and <7 is the Stefan-Boltzmann constant. 

3.3 Description of numerical formulations 

3.3.1 Finite element formulation 

The Galerkin weighted residual method is used to derive finite element formu- 
lations. Within each element, the unknown function T may be approximated at 
any time t by the relationship 

k 

T(x,y,t) = '£N i (z,y),T t (t) (3.13) 

»=1 

where k is the number of nodes assigned to the element, T; are the discrete nodal 
values of T, and Ni are the shape functions. 

To derive the element equations for differential Equation(3.3) and their bound- 
ary conditions, the solution domain R in two dimensions, as shown in Figure 3.3, 
is divided into m linear triangular elements of k nodes each. Application of the 
Galerkin method to Equation(3.3) yields 
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where Ni are the shape functions which are conveniently chosen as weighting func- 
tions in the Galerkin method. 

To reduce the order of the derivatives in the equation above and to introduce 
the influence of the boundary conditions, the second order derivatives are integrated 
by parts 



( 3 . 15 ) 


where i is the unit vector in the X direction, j is the unit vector in the Y direction, 
and n is the outward unit vector normal to the surface, which coincides with the 
boundary to the solution domain. For the surfaces(or sides) of the elements con- 
tained within the solution domain, the surface integrals cancel out by continuity 
of heat flux when the element equations are assembled into the global system of 
equations. Only those surfaces of elements which coincide with the boundary to 
the solution domain, and do not have a specified temperature at the boundary, 
contribute to the surface integral. 

The surface integral can be expressed as the sum of integrals over the external 
surfaces A as follows: 
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( 3 . 16 ) 
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Substitution of Equations(3.15) and (3.16) into Equation(3.14) results in the 
following expression 


dXdY = 


[ QNidS-j h cr (T^-T cr )NidS- ( 0 r (T^ -T r )NidS 
J M J At Ja 4 


(3.17) 


After substituting the expression for unknown function T ^ into Equation 
(3.17), the resulting element equation is expressed as 


f ( dNi dNj dNi dNj\_ ,, , f ( . \ 

U^KexJx + W-9Y) T ‘ iXiY + U ) (<7Wi)«^+ 

[ h„N,N j T j iS+ [ (3 r NiNjTjdS = 

J A& J A 4 

I QNidS+ f h cr NiT cr dS+ [ (3 r NiT r dS 
J a 7 Ja 3 J a 4 


(3.18) 


Finally, the assembly of the nonlinear transient element equations can be ex- 
pressed in matrix form as follows: 


[cm + 


[Kc] + [K h ] + [K r 


{T} = {F q } + {F h } + {F r } (3.19) 


where 


CNiNjdXdY 

(3.19.1) 

„( dNi dNj dNi dNj \ , 


*{ex 3x + ar eY ) dxdY 

(3.19.2) 

■crNiNjdS 

(3.19.3) 


[K r ] = [ PrNiNjdS 

J A* 


(3.19.4) 
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[F q ] = / QNidS 

4 JAt 

(3.19.5) 

[F a ] = [ h cr NiT cr dS 

Ja» 

(3.19.6) 

[F r ] = / PrNiTrdS 

JA t 

(3.19.7) 

3.3.2 Time stepping scheme 


The finite element method is now firmly established to investigate transient 
field problems governed by parabolic equations. However, the time derivative has 
usually been approximated by the finite difference method. Thus, a finite element 
discretization in space is coupled with finite difference approximation in time. 

3. 3.2.1 Imolicit method 


Let t r denote a typical time in the response so that t r+1 = t r + A t where At 
is the time increment, and an intermediate time t$ within each time step may be 
expressed as 

U = t r + 0 < 0 < 1 

(3.20) 

Then, Equation(3.19) at te is given as 


{C]{T} e + mm, = {F} 

(3.21) 

A first-order Taylor expansion in time is used, 


{ry = {r}»- j{t},(»a«) 

(3.22) 


and the following approximation for {T} is introduced. 
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|{r}»= ({ 7 y +1 -{iyW (3.23) 

Substitution of Equation(3.23) into Equation(3.22) yields 

{T} 0 = (1 - 0){T} r +6{T} r+1 (3.24) 

Finally, substitution of the expressions for {T}q and {T}e into Equation(3.21) 
gives 


( a! + {rr+ “ = ( a7 - f 1 " # ) w) < r > r + W < 3 - 25 ) 

Since the values of [C], [K], and {F} depend on {T}$, a choice from among 
the values 6 = 0, 1/2, 2/3, and 1, respectively, yields Euler forward-difference, 
Crank-Nicolson center-difference, Galerkin, and fully implicit backward-difference 
formulations. The fully implicit backward-difference scheme is unconditionally sta- 
ble and predicts a smooth decay. 

3. 3. 2. 2 Explicit methods 

The implicit method is much more stable than the explicit method, but requires 
an iteration within each time step. To avoid iteration, a three level scheme proposed 
by Lees [41] was used by Comini et al. [42-46] and Morgan et al.[47]. Oscillations 
have been observed in certain instances, so a slightly modified form was used to 
improve stability. 

Another three level scheme is referred to as the Dupont three level scheme. 
Hogge[ 48 ] demonstrated its overall performance to be superior in accuracy and 
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stability to other time stepping schemes in solving the one-dimensional homogeneous 
equation. 

Thomas et al.[49] compared several time integration schemes such as the Lees 
scheme, the Dupont scheme, and the Crank-Nicolson method. He concluded that 
the Dupont three level scheme was clearly superior to that of Lees in both ac- 
curacy and stability, and temperature-dependent terms should be evaluated using 
{T} r+1 instead of {T} r+ f. With the use of Dupont method, Equation(3.21) can 
be approximated as follows: 

(jm + ^){rr +! = ^{7T +1 + + m (3.26) 

The Equation(3.26) allows the explicit evaluation of {T} r+2 without iteration, 
provided that {Z’} r+1 and {T'} r are known. However, this scheme is not self- starting, 
so that {T} r+1 at the first time step may be calculated by using the implicit method. 

3.3.2. 3 Latent heat evaluation schemes 

The principal difficulties in the analysis of the phase change problem are that 
the variation of the heat capacity is very severe at the interface, as shown in Figure 
3.2. The position of the moving interface is not known a priori and the shape may 
be multi-dimensional. Thus, physically realistic approximation techniques must be 
used to overcome these difficulties. It is convenient to divide them into two groups, 
based on the choice of grids. 

In the first group, the moving mesh technique continuously tracks the location 
of the interface by deforming the grid system to maintain the finest mesh in the 
vicinity of the critical phase change region. This technique may be limited to very 
simple geometries. In the second group, a fixed grid technique can avoid tracking 
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down the position of the moving interface, but the interface is generally at an un- 
known location between nodes. Many different types of methods are available with 
the fixed grid system. The first method uses the enthalpy as a dependent variable 
along with the temperature, and may be referred to as the enthalpy method[50]. 
Since two dependent variables are used, the system of algebaric equations are solved 
by iteration. 

The second method treats the latent heat effect accompanying a change of 
phase in terms of a temperature-dependent specific heat, or with the use of an 
enthalpy function. These methods avoid the moving interface difficulty so that 
instead of continuously tracking down the position' of the moving interface, the 
same numerical scheme for conduction heat transfer without phase change is equally 
applicable to the phase change region. Then, the position of the interface can be 
easily determined by linear interpolation of nodal temperatures. 

When temperature approaches the phase change temperature, the heat capacity 
tends to the Dirac delta function, and cannot be satisfactorily represented across the 
peak by any smooth function. Frivik and Comini [42] proposed a technique based 
on the integral of the heat capacity with respect to temperature 



This is a smooth function of temperature in the phase change zone. Therefore, 
the enthalpy rather than the heat capacity is interpolated in a element as follows: 

k 

B = (3.28) 

1=1 


where Hi are the enthalpy values at nodal points. 
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From definition, the heat capacity can be expressed as 


C = 


dH 

dT 


(3.29) 


Thus, the values of the heat capacity can be approximated by evaluating the 
gradient of enthalpy with respect to temperature. Defining the direction n to be 
normal to the interface fine, Equation(3.29) is expressed as 


where 


C = 


\dn dn ) 

(<W dH_ \./dr\ 
[dX nx + 8Y ny )'\dn) 


In y 

dr 

dn 


dT dT 

dX / dn 
dT dT 

dY 1 dn 



(3.30) 


Hence, for the entire element, the fined expression of the heat capacity proposed 
by Del-Giudice et al.[51] is given as 


C = 


dH dT 

dxdx 


+ 


dH dT 
dY dY 



2n 


(3.31) 
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CHAPTER IV 

VAPOR FLOW DYNAMICS IN HEAT PIPES 

4.1 Introduction 

Analysis of the hydrodynamics of vapor flow in heat pipes with metallic working 
fluids indicates considerable difficulties at low vapor pressure, due to the extremely 
small vapor densities. Even for relatively small heat transfer rates, vapor velocity in 
the chordwise direction can be very large, accelerated towards sonic velocity because 
viscous action causes pressure and density to decrease. Additional heat input causes 
choking at the condenser inlet. Such behavior makes it necessary to include vapor 
compressibility and viscous action in mathematical models. The vapor pressure 
drop due to friction cannot be recovered completely in the condenser section. Thus, 
the temperature distributions along the length of the heat pipe are not isothermal, 
and thermal resistance in the vapor region is significant. Studies of the distributions 
of temperature, pressure, and velocity in the vapor passage along the length of a 
heat pipe are essential for an evaluation of the maximum heat transfer rates and 
prediction of correct heat pipe performance. 

Reviews of the literature on hydrodynamic processes in vapor flow of cylin- 
drical heat pipes made by Tien[52] and Ivanovskii et al.[ll] indicate that no com- 
pletely detailed investigation has been presented thus far. Moreover, experimental 
measurements of pressure and velocity for metallic working fluids have not been 
reported. 

An analysis of the steady, compressible, one-dimensional, laminar flow of 
sodium vapor is presented for the case of a flat plate-type heat pipe with asymmet- 
rical boundary conditions. In addition, shear stress at the liquid-vapor interface, 
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variations of vapor quality, and momentum and energy factors are considered. A 
similarity solution for a semiporous channel is used to provide the velocity profile 
at cross sections. 


4.2 Similarity solution for a semiporous channel 


4.2.1 Governing differential equations 

A sketch of the geometry for a semiporous channel is given in Figure 4.1. The 
following assumptions axe made: that the fluid flow is incompressible and lami- 
nar, and that the properties of the fluid are constant. The width of the channel 
is assumed much greater than the height. Therefore, two-dimensional flow is con- 
sidered. Fully developed flow is assumed in the channels. A constant injection or 
suction velocity is used. With these assumptions, the Navier-Stokes equations for 
two-dimensional, steady state, incompressible, laminar flow are written as 


„du irm du i dp d 2 u d 2 u 

dX ^ dY pdx ^ [dX 2 dY 2 


(4.1) 


TT av ,,.dv i dp \a l v d u/ -] 

dx dY P dY [ex 2 dY 2 J 


(4.2) 


and the continuity equation is 


dU dV* 
dX + dY 


= 0 


Boundary conditions used on the channel surfaces axe 


(4.3) 


U = 0, 


and V* = 0 at Y = 0 (solid wall) 


(4.4) 
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-Vo : Injection, depicted. 
V e : Suction. 


Figure 4.1. Schematic diagram of a semiporous channel. 
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{7 = 0, and V* = V 0 at Y = D (porous wall) (4.5) 

4,2.2 Similarity transformations 

The governing equations can be transformed into total differential equations 
by the use of a dimensionless length coordinate, u, and a dimensionless function, 
f, which automatically satisfies the continuity equation. The new variables [53] are 
defined as 


u? = — 

D 


(4.6) 


/(«) = 


DV 


(4.7) 


where V(X) = U*(0) - and D r *(0) is the average velocity at X = 0. 

Then, the local velocity can be expressed by the new variables as follows: 


u = w =Vf ' 


(4.8) 


v ° f 


(4.9) 


Use of Equations(4.1), (4.2), and (4.6) to (4.9) gives 


1 dP_ 

pdX 


= V 




(4.10) 


.I2£_v. 

P du> 


Vo ff - pf" 


(4.11) 
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Since the right side of Equation(4.11) is a function of u; only, differentiating 
with respect to X yields 


8 2 P 

dXdu 


(4.12) 


On taking the derivative with respect to u >, applying Equation(4.12), and inte- 
grating, Equation(4.10) becomes 


Re 0 


(/') 2 - //" 


+ /"' = -A 


where Re 0 = and A = — — 14 

V l/p O J\ 

The transformed boundary conditions axe written as 


(4.13) 


f' = f = 0 at a? = 0 (solid wall) 


(4.14) 


f' = 0, / = 1 at u> = 1 (porous wall) (4.15) 

4.2.3 Solution of equation 

Differential Equation(4.13) together with the associated boundary conditions, 
constitutes a nonlinear boundary value problem with the parameter Re 0 ■ Since this 
governing equation is no longer partial differential equation, it can be solved numer- 
ically by the Runge-Kutta integration method. However, the number of boundary 
conditions is not enough to solve Equation(4.13). Therefore, for each specified Re 0 , 
the value of A and /"( 0) are guessed to solve a system of the first order differen- 
tial equations. Then, at u> = l(porous wall), the calculated values of f and / are 



48 


compared with the given boundary conditions. When both values are not accept- 
able, other values are used until the convergent tolerance is satisfied. When wall 
Reynolds numbers for suction are greater than 13, separation occurs on the solid 
wall. Similarity solutions are valid for wall Reynolds numbers up to 13. Hence, the 
equation is solved for wall Reynolds numbers ranging from - 30 to 13. 

A comparison of the present results with those found in the literature[53,54], 
which use the perturbation method for small Reynolds numbers and a different 
numerical scheme, verifies the accuracy of this numerical solution. 

The velocity profiles for the suction and injection sections are shown in Figures 
4.2 and 4.3. Asymmetric boundary conditions cause the velocity profiles to be 
asymmetric. For injection, Re 0 < 0, the location of the maximum velocity shifts 
from the center of the channel toward the solid wall. Thus, injection increases 
friction at the solid wall, and decreases friction at the porous wall. As wall Reynolds 
numbers increase, the degree of shifting is large, although the general shape is 
changed little. For suction, Re 0 > 0, the location of the peak velocity shifts toward 
the porous wall. Hence, suction increases friction at the porous wall and decreases 
friction at the solid wall. Unlike injection, separation for suction appears around 
a wall Reynolds number of 13. The velocity profiles change considerably with wall 
Reynolds number. In general, friction for the semiporous channel is larger than that 
for the impermeable channel. 

Results[19] from two-dimensional, incompressible, Navier-Stokes equations for 
a flat plate heat pipe show similar velocity profiles. In the evaporator, the velocity 
profiles retain similarity for a wall Reynolds number of 50, which is the highest 
value tested. However, back flow was observed near the end of the condenser for 
a wall Reynolds number of 10. For a wall Reynolds number of 50, back flow was 
observed in the entire condenser section. 
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4.3 Compressible vapor flow analysis 

Steady, compressible, one-dimensional, laminar flow in a heat pipe is consid- 
ered. The principal governing equations for mass, momentum, and energy axe 
formulated by using the average velocity as described in Chapter 2. This veloc- 
ity is approximated from velocity distributions based on the similarity solution of 
semiporous channels. Shear stress at the interface, and the momentum and energy 
factors are similarly calculated and shown in Figure 4.4. 

The fluid in the vapor passage of the heat pipe is assumed to be a mixture 
of liquid and monatomic vapor. Thus, the quality of the vapor is considered. The 
specific volume v and the enthalpy h axe expressed as 

V = Vf + X q x (Vg - Vf) (4-16) 


h — hf -f* Xq X hfg 


(4.17) 


The specific volume of the saturated vapor can be approximated by 


R U T 

Vg ~ PM 


(4.18) 


Also, the temperature and pressure are related by the Clausius- Clapyron equa- 
tion: 


dP h fg M dT 

P R u T 2 


(4.19) 


where hf g is the enthalpy of vaporization, Vf and v g are the specific volume of 
saturated liquid and vapor, respectively, and X g is the vapor quality. 

ORIGINAL PAGE IS 
OF POOR QUALITY 



ORIGINAL PAGE IS 
OF POOR QUALITY 


lo^om iCSjaua JO uiTV}iiamofl 

M 


01 o 
— 4 -iri 



w Sh 
0 ) 

•§ 
9 6 
>«n d 

1 £ 


OT 

. B 

'd o 

o o 
cocu 


w 

0 *d 

21 

1 CJ 

0 ) 

, 5 « 




<d <d 3 
§§■ 8 * 

o o » a) 

££££ 

n n a n 


• o □ ■ 


o> oz 

ips^i pn°S Y 8 Jo^o^i uoipT-y 


o-n- o st- 0 6I- 0 G2- 0 i2- o ie- 
IT 8 ^ snoaoj ys jo^o^ xioipi-y 





Figure 


53 


4.3.1 Formulation of differential equations 

Bankston and Smith[17] showed that variation of the momentum and energy 
factors with axial distance and radial Reynolds number is very small except near 
the end of the heat pipe. Figure 4.4 shows that both Mf and Ef are nearly indepen- 
dent of the radial Reynolds number except near separation. These results are taken 
from the similarity solution for semiporous channels. Therefore, it is assumed that 
the derivatives of Mf and Ef with respect to chordwise distance axe equal to zero. 
In addition, assuming Vf and hf g are constant quantities, and combining Equa- 
tions^. 15) through (2.21) and Equations(4.16) through (4.19) yields the chordwise 
gradients for the density, quality, velocity, pressure, and temperature. 

4. 3. 1.1 Density 

The differentiation of Equation(4.18) with respect to S gives 


dv g _ Vg dP Vg dT 

ls-~PlS + TdS 


(4.20) 


Substitution of Equation(4.19) into the equation above yields 


dVg _ Vg 

dS P 


RuT 

hfgM 


- 1 


dP 

dS 


(4.21) 


The derivative of the mixture specific volume in the S-direction is written as 


dv dv f x dX q 

dS ~ dS + ~ dS q 


dVg 

dS 


dvf 

Is 


(4.22) 


With the prescribed assumption for v f , substitution of Equation(4.21) into 


Equation(4.22) results in 
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dv 

dS 


, \ dX q . v g^-q 

( V °- V f ) ~df + J p- 


RuT 

hfgM 


~ 1 


dP 

dS 


(4.23) 


The derivative of specific volume and density are related as follows: 


= (4.24) 

dS v 2 dS 

From Equations(4.23) and (4.24), the expression for density in differential form 
can be written as 


dp _ 1 

dS = ~v* 

4.3.1 .2 Quality 

From Equation(4.23), the quality gradient is expressed as 


^dXg , VgXq ( R U T ,\dP 


(4.25) 


dX q _ _ 

1 

' dv VgXg ( 

RuT > 

UP' 

dS 

(vg ~ v f ) 

ds + p V 

hfgM, 

1 dS 


Substitution of Equation(4.24) and the equations for conservation of mass, 
Equation(2.15), and momentum, Equation(2.16), into Equation(4.26) yields 


dXq VZ f [ 1 | ***« ( 

dS - {v g -Vf)V 2 \ [ Mf P v 2 \ 
v 2 2 m 0 


RuT_\]dP 

hfgMJl dS 


■FV 2 \ 

SDVM '> (4.27) 


(Vg - Vf) DV 
4. 3. 1.3 Velocity 

The equation for conservation of mass,Equation(2.15), can be expressed as 


follows: 
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dV vrh 0 V dv 

1s = ~d~ + 


(4.28) 


Substitution of Equation(4.23) into the expression above gives 


dV vm 0 V(v g - v t ) dX q VX q v g / R U T \ dP 
dS ~ D + v dS + P v \hf g M 7 dS 


(4.29) 


4.3. 1.4 Pressure 

With the previously mentioned assumptions, substitution of the energy Equa- 
tion^. 17) into the momentum Equation(2.16) yields 


dP Mfpdh Mfrh 0 /, E f V 2 F 0 2 \ FpV 2 

dS E f dS EfVD 7 2 2 ) ~ 8D 

In order to obtain an expression for the derivative of the vapor enthalpy, dif- 
ferentiation of Equation(4.17) is taken with respect to S, and Equation(4.19) is 
substituted: 



dh dX q c p R u T 2 dP 
dS ~ u dS + hf g M P dS 


(4.31) 


After the expressions for the derivative of the enthalpy and quality are substi- 
tuted into Equation(4.29), the pressure gradient can be expressed as follows: 


dP 

dS 


Mf m 0 

~ E f VD 


2hfg + 




J_ hta-Jl v a~ v t EXl 

E t v 8 D v 2 SD 


+ 


Ef V 2 


Ef 




1 - 


h /g M 


M f v a -Vf CjMu 

Ef V 1 hfgM P 


v 


( 4 . 32 ) 


4.3. 1.4 Temperature 


The temperature gradient is derived from the Clausius-Clapyron equation, 


Equation(4.19), as follows: 
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CHAPTER V 

COMPUTATIONAL PROCEDURES 



The algorithm proceeds similarly to most finite element method procedures. 
First, the element data are generated by a grid generation program. This data 
consists of the coordinates cf the nodal points of each element, by element, and 
the properities and boundary conditions for each element. Based on this data, the 
capacitance matrices, [C], and the conductance matrices, [K c ], related to the time 
derivative of nodal temperatures and conduction, are calculated for each element. 
The conductance matrices, [Kh] and [K r \, and the vectors {.F g }, {i 7 /,}, and {F r } 
for specified surface heating, convection, and radiation boundary conditions, are 
computed only for elements having the boundary conditions above. These boundary 
conditions may be time dependent, so these matrices and vectors must be evaluated 
based on the appropriate temperatures. For an implicit time step, the iteration 
scheme is such that temperatures axe initially assumed, and the proper value of 6 
is chosen to obtain stable and accurate results. In this study, 0 = | or 1 is used, 
so iterations are required within each time step. However, for an explicit time step, 
such as the Dupont scheme, the element temperatures for the previous step are 
involved in obtaining the matrices and vectors related to the boundary conditions. 
Then, these element matrices and vectors are assembled into the global matrices for 
the entire solution domain, and the full set of equations may be written in matrix 


form as 
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AM 


{T} = {RM} 


(5.1) 


where [AM] represents the banded coefficient matrix, {T} denotes a column ma- 
trix of unknown temperatures, and {RM} is the column matrix of constants. If 
temperatures are specified on some boundary surfaces, these boundary conditions 
must be incorporated into the global matrices as described in reference [55]. 

Once the global matrices have been assembled, they are solved by using Chole- 
ski decomposition. Choleski’s method[56] has the advantage of being simpler and 
easier to implement than other elimination methods. It may also be used to ad- 
vantage for storage in the computer by overlaying the upper and lower triangular 
matrices in the same storage location([AM] matrix). 

The system of Equations(3.25) for the implicit method or (3.26) for the explicit 
method is then marched forward in time. Among the various iteration schemes, 
the Newton-Raphson method is used for the implicit method because of its fast 
convergence. For the explicit method, iteration is not required, but it is not a self- 
starting scheme, so the implicit method is employed for the first few time steps. 
When temperatures at two previous time steps are known, the explicit method can 
be marched forward. The algorithm flowchart is shown in Figure 5.1. a. 

5.2 Compressible vapor flow 

The five dependent variables: density, quality, velocity, pressure, and tempera- 
ture, are coupled by differential Equations(4.25), (4.27), (4.29), (4.32), and (4.33), 
which are first order, nonlinear, ordinary differential equations. A computer code 
has been written to numerically solve these equations simultaneously using the 
Runge-Kutta integral method, which needs only proper boundary conditions for 
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the first step of integration. 

Since the velocity term appears in the denominator of Equation(4.32), the phys- 
ical boundary condition for velocity cannot be used directly. To avoid this problem, 
proper boundary conditions at the upstream end of the evaporator are determined 
for initiating all calculations as follows: a new boundary condition for the velocity 
is determined a short distance away from the beginning of the evaporator, assum- 
ing incompressible and saturated vapor flow at the temperature corresponding to 
the heat pipe operating temperature. Saturation pressure is used for the pressure. 
Since the velocity is small at this point close to the upstream end of the evaporator, 
boundary conditions determined this way axe realistic. The differential equations 
are solved using these boundary conditions. 

5.3 Coupling vapor flow effects 

When the Mach number is less than about 0.2, compressibility is neglected, 
and friction effects at the interface may be negligible due to the low velocity of 
the vapor. Thus, most studies of heat pipe performance [5, 23-25] assumed that the 
temperature is uniform throughout the vapor space. This approximation eliminates 
difficulties encountered in solving the vapor flow dynamics, and gives simple results 
for low operating temperatures with small heat fluxes. When the vapor temperature 
is uniform in the vapor space, no thermal resistance exists, so a maximum amount of 
energy can be transferred through the vapor space for a given operating condition. 

The energy stored in the vapor space is negligible due to low density, so the 
energy entering the evaporator is equal to the that leaving the condenser. Equa- 
tion(2.12), which describes evaporation and condensation, is applied to every ele- 
ment at the interface, and the energy balance in the vapor space is: 
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E 




mc r p* 



(5.2) 


where m e is the number of elements at the interface in the evaporator, m c is the 
number of elements in the condenser, and A £, is the side length of an element at 
the interface. Equation(5.2) can be solved by iteration to obtain the uniform vapor 
temperature, as long as the temperature at the interface is known. The heat flux 
at the interface is evaluated using Equation(2.12) with known temperatures. 

However, when liquid metal is used as the working fluid, and the heat flux 
applied to the surface is large, the velocity of the vapor is large, so the temperature 
drop in the vapor space should be considered. The temperature drop in the vapor 
region implies that the thermal resistance allows less energy to be transferred from 
the evaporator to the condenser. Thus, the evaporator temperature is higher, and 
the condenser temperature is lower, than in the uniform vapor temperature case. 

When vapor flow dynamics is coupled with the heat pipe shell and capillary 
structure at the interface, the governing equations are solved simultaneously with 
unknown boundary conditions at the interface. Iterations axe required for every 
time step until both results match at the interface. This method may yield accurate 
results, but in general consumes much computational time, due to iteration. Also 
convergence may not be reached. Hence, an approximation method is employed to 
eliminate these difficulties. 

Instead of simultaneously solving governing equations for both regions, govern- 
ing equations for the vapor region are separately solved with a given heat flux, so 
that the temperature drop,A!F 3 , can be obtained. The thermal resist ance,i? ? , in 
the vapor space may then be evaluated from Equation(5.3): 
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Q = 


AT, 

Rg 


(5.3) 


In order to obtain the thermal resistance, the total heat flow rate is needed. 
But without solving both governing equations, the heat flow rate is also unknown. 
Thus, the known heat flow rate at the previous time step is used, and the heat flow 
rate evaluated at the present time step is used at the next time step. 

Coupling of the evaporator and condenser sections becomes very difficult when 
vapor temperature changes along the length of the vapor section. An approximate 
coupling method is used here which accounts for thermal resistance in the vapor 
region but which uses a constant vapor temperature. Since the total heat transfer 
rate through the vapor space, which is computed by using Equation(2.12), is the 
same as the heat input at the liquid-vapor interface in the evaporator or heat output 
at the interface in condenser, a certain fictitious layer, which has the same thermal 
resistance as that evaluated in the vapor space, may be placed at the liquid-vapor 
interface. Thus, the temperature at the new interface in the evaporator is lower 
than that at the true interface, due to the resistance in the fictitious layer, and 
vice-versa for the condenser. With this coupling technique the heat fluxes are 
correct throughout the pipe and vapor thermal resistance is taken into account. 

In order to compute the uniform vapor temperature, the temperatures at the 
new interface are used for Tf { in Equation(5.2), and the heat flux at the interface 
is calculated by using Equation(2.12) with the uniform vapor temperature and the 
temperatures at the new interfaces. This heat flux then serves as the boundary 
condition at the liquid-vapor interface when solving the governing equations for the 
heat pipe shell, capillary structure, and vapor flow at the next time step. By using 
this approximation scheme, iterations can be excluded, while vapor flow effects are 
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retained. 


5.4 Overall computational procedures 

Since the working fluid of the heat pipe is initially in the solid state and the wick 
structure is saturated, the transient conduction equation is applied to the heat pipe 
shell and wick structure. The variable heat flux and radiation boundary condition 
are considered on the outside surface of the heat pipe. The adiabatic boundary 
condition is used at the liquid-vapor interface, due to the vaccuum in the vapor 
space, until the nodal temperatures at this interface are greater than the transition 
temperature (700 K). The implicit method is used for the the first ten time steps, 
and then the explicit method is employed using a time step of 10 seconds. 

When the temperatures of the first two nodes at the liquid-vapor interface are 
greater than 700 K, Equation(5.2) is solved for the vapor temperature by using 
the Newton- Raphson method. Then, the heat fluxes are calculated using Equa- 
tion(2.12) with the known interface and vapor temperatures. These heat fluxes are 
to be used as boundary conditions at the interface for the next time step. However, 
the adiabatic boundary condition is still applicable for the rest of the interface. 
From then on, a small time step is used to obtain stable startup. These procedures 
continue until the heat transport in the chordwise direction is less than the sonic 
limit. 

When the sonic limit is encountered, the total heat transport through the 
vapor space should equal the sonic limit, and Equation(2.14) is solved for the vapor 
temperature. The heat fluxes at the interface in the evaporator are calculated by 
using this vapor temperature and the nodal temperatures at the interface. Then 
the heat fluxes for the condenser are evaluated in proportion to nodal temperatures 
which are greater than 700 K. Again, these heat fluxes are used to solve the transient 
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conduction equation for the heat pipe shell and wick structure. This scheme is 
effective until the entire vapor space achieves continuum flow and the heat transport 
in the chordwise direction is less than the sonic limit. 

When the sonic limit is not encountered, Equations(4.25), (4.29), (4.32), and 
(4.33), which are first order, nonlinear, ordinary differential equations for compress- 
ible vapor flow dynamics, are solved to evaluate the thermal resistance in the vapor 
space, as stated in section 5.3. This computational scheme is used until the desired 
operating temperature, or the steady state condition, is reached. Figures 5.1. a and 
5.1.b schematically shown the computational procedures. 

Since any implicit scheme is excluded, except the first few steps to provide 
initial conditions for the explicit method, time steps are gradually decreased from 
0.5 seconds to the order of milli-seconds to maintain a stable condition at the liquid- 
vapor interface. The computer code is successively tested against the given physical 
model by following these procedures. 
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CHAPTER VI 

RESULTS AND DISCUSSION 

6.1 Transient conduction and phase change 

To verify the finite element formulation and algorithm, computed numerical 
solutions are compared with available solutions, such as analytical and other ap- 
proximate solutions. 

6.1.1 Temperature of a semi-infinite body 

This case involves pure conduction of heat, without phase change, over a semi- 
infinite body. It is solved as a two-dimensional problem, shown in Figure 6.1. 
Adiabatic boundary conditions are assumed throughout, but at the surface (X = 
0), the constant heat flux ( Q = 100 kW/m 2 ) is imposed to heat the surface while 
the initial sodium temperature of 293 K is uniform. The conductivity and specific 
heat are assigned constant values for sodium. 

The exact solution to this problem was given by Luikov[57] as 


T(X,t) = ?g-Vrt 



2Vai er/C (2v / Si). 


T, 


( 6 . 1 ) 


where a = — 
pc P 

For the numerical solution, an explicit time stepping scheme is used with an 
implicit method for starting. A time step of 5 seconds is used, and numerical 
calculations are terminated as the temperature of the last node starts to change 
from the initial temperature. Numerical results are compared to the analytical 
solution in Figure 6.2. The results yield excellent agreement. 
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Figure 6.1. Finite element mesh for the transient conduction heat transfer problem. 
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6.1.2 Convection and radiation boundary conditions 

This example tests the ability of the numerical method to handle convection 
and radiation boundary conditions. The same grid system shown in Figure 6.1 is 
used, but the boundary condition on the surface at X = 0 is replaced. For a Biot 
number of 1.0, three different cases axe examined. First, only a convection boundary 
condition (/3 = 0.0) is applied to the surface, and an initial temperature of 1650 K 
is used. Next, convection and radiation boundary conditions(/3 = 1.0) are applied 
to the surface with the same initial temperature. Finally, an initial temperature 
of 2650 K, which is greater than the previous two cases, is used for (3 = 4.0. In 
these tests, phase change is not involved, and the properties of sodium are used for 
specific heat and conductivity. 

An analytical solution for this case is not available except for one having only 
a convection boundary condition, so the numerical results are compared with prob- 
ablity solutions which Sheikh and Sparrow[58] found for the same cases. As shown 
in Figure 6.3, temperatures at the boundary surface and opposite end are plotted 
versus time. The results show good agreement. 

6.1.3 Phase change of sodium 

In this example, sodium is used as the phase change material for two-dimens- 
ional region as shown in Figure 6.4. An initial temperature of 393 K is used. The 
ratio of diffusivities of solid and liquid phases is assumed to be one to compare 
with the existing solution. Boundary condi tions(T 4 = 293 K) of the first kind 
are imposed on the surfaces. Since specified temperature boundary conditions are 
suddenly imposed on the surface, the temperature gradient in this region is infinite, 
which causes unstable conditions. This difficulty may be eliminated by introducing 






FICTITIOUS LAYER 


ORIGINAL PAGE IS 
OF POOR QUALITY 


71 


ADIABATIC CONDITION 

0.3 U 



Figure 6.4. Two-dimensional mesh used to represent a corner region; 242 elements. 


144 nodes. 


ADIABATIC CONDITION 




72 


a layer of fictitious elements with negligible thermal storage capacity and very high 
conductivity, such that the temperature gradient at this region artificially becomes 
finite. For this purpose, numerical values of the thermophysical properties of the 
ficititious layer are chosen as follows: 

pc p = 1.26 x 10 4 J/m z K 
K = 4.19 x 10 3 W/mK 

Two different time stepping schemes, the fully implicit method and the explicit 
method, are employed to compare the computational time and results. A time step 
of 5 seconds is used. For the two-dimensional case, the position of the interface 
is compared with an approximate solution given by Rathjen and Jiji[59] . This 
solution assumes that the ratio of diffusivities of solid and liquid is unity, and 
that the interface at points beyond three times of the one-dimensional interface in 
X and Y directions is the same as the one-dimensional interface position. Both 
numerical solutions for explicit and implicit methods are directly compared with 
the approximate solution, and good agreement is obtained for early time steps, as 
shown in Figure 6.5. As time passes, the position of the interface given by numerical 
results advances further than for the approximation, and the interface lines obtained 
from the approximate are not perpendicular to the boundary surface. According 
to assumptions made by Rathjen and Jiji, adiabatic boundary conditions for the 
approximate solution can be located beyond the present solution domain at certain 
time steps. However, the numerical calculations use adiabatic boundary conditions, 
so that all heat is used for phase change, while the approximate solution does not 
match the same boundary condition within the solution domain. This explains why 
the positions of the two interfaces are different. 

Numerical results predicted by the explicit method using an implicit scheme for 
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a few initial time steps are close to those of the fully implicit scheme and the explicit 
method consumes much less computational time, since iterations are eliminated. 

Even though the numerical solutions do not exactly match the approximate 
solutions, because of different boundary conditions, general trends in the numerical 
results and the usefulness of the explicit time stepping scheme are verified. 

6.2 Compressible vapor flow 

To test the differential formulations and algorithm for vapor flow, the numerical 
solutions are compared to other numerical results[18] for cylindrical heat pipes with 
sodium as the working fluid, since data for a flat plate heat pipe axe not available. 
The two-dimensional mass, momentum, and energy conservation equations, which 
were transformed in terms of the stream function, vorticity and enthalpy, were 
solved numerically. The dimensions of this cylindrical heat pipe are: total length 
of 0.6 m, evaporator length of 0.2 m, condenser length of 0.3 m, and inside radius 
of 0.0086 m. Since variation of temperature and pressure in the axial direction for 
low heat input is small, two different operating conditions are selected. One has 
an operating temperature of 818 K with a uniformly distributed heat input of 610 
watts. Another has an operating temperature of 841 K with a heat input of 1265 
watts. 

To match the results for the rectangular heat pipe to that for the cylindrical 
heat pipe, the rectangular cross-section of the vapor space is modeled such that the 
cross-section area is the same as for the cylinder, and the hydraulic radius is close 
to the radius of cylinder. Also, the same heat input is applied to yield the same 
velocity and a similar Reynolds number in the axial direction, and the same length 
for each section is used. The wall Reynolds number at the condenser must be less 
than 14 to use the results of the similarity solution. Therefore, two different cross 
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sections are chosen. For high heat flux, the wall Reynolds number cannot be less 
than 14 so that a square of 0.0152 m is used and comparison is made only for the 
evaporator and adiabatic sections. For low heat flux, the cross section is chosen such 
that the wall Reynolds number at the condenser is close to the allowable maximum 
value. Thus, the cross section is rectangular with a width of 2.11 cm and a height of 
1.1 cm, resulting in a wall Reynolds number of 13.4 in the condenser. Even though 
the two models do not match exactly, the effect of vapor flow may be simulated 
using these models. 

One- dimensional numerical results of the flat plate heat pipe are compared 
with published two-dimensional numerical data for cylindrical heat pipes as shown 
in Figures 6.6 and 6.7. Pressure and temperature drops for the one-dimensional case 
are a little greater than those for the two-dimensional case, but overall agreement is 
good. For high heat flux, pressure variation with axial position is somewhat differ- 
ent for the two cases, but the total pressure drops are nearly the same. Temperature 
profiles agree well. Thus, as shown in Figures 6.6 and 6.7, one-dimensional differen- 
tial formulations yield valuable results. The difficulties of solving two-dimensional 
governing equations are avoided. 

To investigate various vapor effects, two other heat pipes with sodium as the 
working fluid were selected for analysis. One has an adiabatic section, the other 
does not. Dimensions of the heat pipes are shown in Table 6.1. Different operating 
temperatures were chosen so that temperature effects could be evaluated. For each 
operating temperature, various heat fluxes were considered. Since the similarity 
solutions used are valid only while wall Reynolds numbers are less than 14, the 
maximum heat flux on the condenser is chosen to satisfy this condition. 
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Figure 6.6. Comparison of axial variation of vapor pressures. 
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Figure 6.7. Comparison of axial variation of vapor temperatures. 
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Figure 6.8 shows solutions for axial variation of vapor temperature, pressure, 
velocity, and density, as obtained from Equations(4.25), (4.27), (4.29), (4.32), and 
(4.33). An operating temperature of 773 K and a uniformly distributed heat transfer 
rate of 200 watts were used. Momentum equation(2.16) implies that variation of 
pressure in the axial direction depends on the relation between the contributions of 
inertia and friction. Thus, the pressure of the vapor in the evaporator falls sharply 
along the vapor passage due to friction and acceleration of the flow caused by the 
injection of mass. The corresponding temperature also drops sharply about 5.5 K. 
Since the density of sodium vapor is relatively low at low temperatures, the velocity 
is correspondingly large in order to transfer the required mass. Results show that 
the maximum Mach number reaches 0.3 at the exit of the evaporator. 

In the condenser, extraction of mass tends to increase pressure because of 
decreasing velocity while friction tends to decrease pressure. Therefore, pressure 
cannot recover completely due to friction loss. Even though some temperature 

m 

recovery is achieved in the condenser, the temperature at the end of the condenser 
is about 4.2 K less than that at the beginning of the evaporator. Figures 6.9, 6.10, 
and 6.11 show variations of temperature, pressure, and Mach number, respectively, 
corresponding to five different heat fluxes. At an operating temperature of 773 K, 
larger heat input leads to greater pressure and temperature drops, and to a higher 
Mach number. When the Mach number exceeds about 0.2, variation of the Mach 
number is large, due to expansion of the vapor. For heat inputs of 50 and 100 watts, 
the vapor can be assumed to be isothermal. 

Figure 6.12 shows that heat transfer is limited to about 351 watts, because 
velocity at the end of the evaporator approaches the sonic velocity at this heat 
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Figure 6.9. Axial variation of vapor temperature for various heat inputs. 






ORIGINAL PAGE IS 
OF POOR QUALITY 


82 


I 

I 

I 



Figure 6.10. Axial variation of vapor pressure 
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at the sonic limit condition for operation temperature of 773 K. 
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Figure 6.13 presents axial variation of vapor temperature, pressure, velocity, 
and density at an operating temperature of 808 K, with a uniformly distributed 
heat input rate of 350 watts. For this case, the operating temperature is 35 K 
higher, but the pressure and density are twice as high as for 773 K. Thus, velocities 
are small and the maximum Mach number is only about 0.2. This leads to a small 
pressure drop, with a corresponding temperature drop of about 2 K. Even though 
the heat input is twice as much as for 773 K, a comparison of temperature drops 
suggests that the pressure and temperature drops depend mainly on the operating 
temperature. 

In the adiabatic section, there is no mass injection or extraction. Hence, pres- 
sure simply decreases due to friction, and the velocity increases very little. The 
contribution of the adiabatic section is the addition of pressure and temperature 
drops such that, for a high Mach number at the exit of the evaporator, maximum 
heat transfer rate is reduced considerably. Figure 6.14 shows this phenomenon. A 
heat input of 865 watts results in a Mach number of 0.7 at the exit of the evapora- 
tor. Then, the vapor is accelerated in the adiabatic section, so that velocity of the 
vapor approaches sonic velocity at the end of this section. 

As expected, in the condenser section the pressure and temperature are par- 
tially recovered, but the degree of recovery is small due to the small vapor velocity 
and large friction loss in the longer condenser. Numerical results show that temper- 
ature recovery is only 0.6 K. Variations of temperature, pressure, and Mach number 
corresponding to six different heat fluxes at an operating temperature of 808 K are 
presented in Figures 6.15, 6.16, and 6.17, respectively. For heat input up to 200 
watts, pressure recovery is not observed. This implies that the pressure drop due to 
friction in the condenser may exceed the value of the inertial contribution. However, 
total pressure drops are small, so that the vapor can be assumed isothermal. 
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at the sonic limit condition for an operating temperature of 808 K. 







Figure 6. 15. Axial variation of vapor temperature for various heat inputs. 









Figure 6.16. Axial variation of vapor pressure for various heat inputs. 








Figure 6.17. Axial variation of Mach number for various heat inputs. 
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Figure 6.18 shows the distribution of nonuniform heat flux on the evaporator 
length of 8 cm and condenser in another case studied numerically. Axial variations of 
vapor temperature , pressure, velocity and density corresponding to the distribution 
shown in Figure 6.18. a, are presented in Figure 6.19. At a short distance from 
the beginning of the condenser, the slight heat extraction causes that the axial 
variation of the vapor velocity is small. Thus, the friction effect in this region 
is dominant. The minimum temperature and pressure appear in the condenser, 
instead of at the exit of the evaporator. As heat extraction increases, the absolute 
value of axial variation of velocity increases and the vapor velocity decreases. The 
pressure is recovered gradually. In Figure 6.20, corresponding to the distribution 
shown in Figure 6.18.b, temperature and pressure are recovered immediately from 
the beginning of the condenser, due to large heat extraction. As expected, axial 
variations of temperature, pressure, velocity, and density depend on distribution of 
heat input. 

In summary, the degree of pressure recovery depends mainly on the heat flux, 
operating temperature, and length of the condenser. The higher the heat flux, the 
lower the operating temperature and the shorter the condenser, the greater the 
pressure recovery. However, for this case, the pressure drop in the evaporator can 
be large enough so that the heat pipe is not entirely isothermal. Since experimental 
data is not available, quantitative comparison cannot be achieved. However, com- 
parison of the general behavior of the present results with published data[15,18] 
gives qualitative agreement. 
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Figure 6.13. Distribution of heat input to the evaporator and condenser. 
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Figure 6.19. Axial variation of vapor temperature, pressure, velocity, and density 
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Figure 6.20. Axial variation of vapor temperature, pressure, velocity, and density 
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6.3 Heat pipe startup 

To verify the model developed for heat pipe cooled leading edges, computed 
numerical results are compared to Camarda’s experimental results[30]. The shape 
of the cross-section of Camarda’s model is different from that of the model un- 
der consideration. However, the rectangular cross-section used in this research is 
modeled such that the cross-section areas occupied by the heat pipe shell, the cap- 
illary structure, and the vapor space are the same as for Camarda’s experiment. A 
comparison of these two cross-sections is shown in Figure 6.21. The surface area 
where the heat flux is applied for the model is chosen to be identical to that of 
the experiment so that the width and length of the heat pipe are identical, but the 
thicknesses of the components of the two heat pipes are different. The thickness is 
much smaller them other dimensions, and the temperature drop at the cross-section 
is small compared to that in the axial direction, so the effects of the difference in 
thicknesses are minimal. Effects of energy storage and fusion of the working fluid in 
the heat pipe shell and the capillary structure are the same, and the effect of vapor 
flow dynamics may be simulated. Since the minimum dimension of the vapor space 
is different, due to the different thickness of the heat pipe, the transition tempera- 
ture, which is evaluated based on the dimension of the vapor space, is not identical. 
However, this temperature determines the boundary condition at the liquid-vapor 
interface. Therefore, to simulate Camarda’s case, the inner diameter of the circular 
heat pipe is used to calculate the transition temperature. 

The material of each component used in the computation is the same as in the 
experiment, except that the braze alloy is assumed to be Hastelloy X to simplify 
computations in the heat pipe shell. Properties of the materials are varied with 
temperature during numerical calculation. 
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Figure 6.21. Comparison of the two cross-sections. 
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Heat distribution on the heat pipe surface as used in the model is based on a 
normalized distribution[30j is shown in Figure 6.22 and heat input at the stagnation 
point is shown in Figure 6.23. Actual heat distribution due to aerodynamic heating 
is shown in Figue 6.24. In addition to aerodynamic heating, a radiation boundary 
condition is used to extract energy from the entire heat pipe surface, so that the net 
rate affects the heat pipe internal operation. The section which has a net positive 
rate is considered the evaporator, and the section with the net negative rate is the 
condenser. 

Figure 6.25 shows the two-dimensional grid system, dimensions, boundary con- 
ditions, and materials used to represent a leading edge. It is assumed that the 
temperature of the heat pipe is initially the ambient temperature, which is below 
the melting point of the sodium working fluid. A temperature of 700 K is used for 
transition from free molecular to continuum flow. 

An explicit method with a few implicit initial steps is used as the time stepping 
scheme. Until continuum flow is established in the vapor space, a constant time 
step of 10 seconds is used. This time is gradually decreased, due to the increasing 
heat input on the external surface. 

Figure 6.26 shows a comparison of numerical and experimental results at var- 
ious times. Since measured temperatures were not surface temperatures but the 
temperatures at intermediate points, numerical data at the interface between the 
heat pipe shell and capillary structure are used for comparison. Since most of the 
energy which is transferred by conduction through the heat pipe shell and capillary 
structure in the axial direction is added near the stagnation point, temperature at 
that point increases rapidly and a large temperature gradient is established in the 
heat pipe during startup. This implies that the heat pipe is not effective, due to 
the extremely small vapor density, for about the first 450 seconds. Temperatures in 







Figure 6.23. Heat input at the stagnation point. 




Heat Flux Distribution, kw/m 2 

0.0 50.0 100.0 150.0 200.0 250.0 



Length of Heat Pipe, s / l 


Figure 6.24. Heat flux distribution on a leading edge model 
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Figure 6.26. Temperature of a heat pipe cooled leading edge 
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the rest of the heat pipe are not much different from the initial condition. At 550 
seconds, temperature at the stagnation point is greater than the transition temper- 
ature and the temperature gradient in the region adjacent to this point begins to 
decrease. This means that a continuum flow region is established in the vapor space, 
so that energy is mainly transported via latent heat transport. However, as shown 
in Figures 6.27 and 28, most of the working fluid is in the frozen state, and free 
molecule flow prevails in most of the vapor space. With time, the continuum flow 
region is expanded, so that at about 750 seconds, two-fifths of the device performs 
as a heat pipe. The rest of the heat pipe still behaves like a solid bar in which heat 
is transported only by conduction, even though the working fluid is melted over the 
entire length. Between the two flow regions, a large temperature gradient exists, 
where the continuum flow front is located. This front moves down toward the cool 
end of the heat pipe. At 1000 seconds, most of the heat pipe is active,, except near 
the end of the condenser section. About 100 seconds later, continuum flow exists in 
the entire vapor space as shown in Figure 6.28. Until this moment, the sonic limit 
is encountered due to the large temperature gradient in the vapor space. Therefore, 
maximum heat transport in the vapor space is equal to the sonic limit. However, 
the expression used for the sonic limit excludes friction effects at the interface, so 
that the sonic limit is overestimated. Thus, at a time of 1000 seconds, temperature 
near the stagnation point is less than that found in experiments, and the continuum 
flow front is advanced further. After a time of 1200 seconds, the heat pipe becomes 
isothermal, but due to increasing heat input the temperature increases until the 
system reaches a steady state at a later time. 

Figure 6.29 shows temperature history of the heat pipe at three different loca- 
tions in the axial direction. When a part of the vapor space has free molecular flow, 
most of the heat input is used to heat up the condenser section, so that the temper- 
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Figure 6.28. Position of continuum flow front. 
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Figure 6.29. Temperature history of three locations in the heat pipe. 
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ature increases very slowly. As the continuum flow front approaches the end of the 
heat pipe, the remainder of the vapor space quickly reaches continuum flow. After 
the entire heat pipe is active, temperatures at all axial locations increase relatively 
fast, and then tend to change slowly as steady state is reached. Experimental results 
exhibited steady state at a time of 1500 seconds where the maximum temperature 
was 883 K. However, the numerical results yield a maximum temperature of 948 K 
at 1450 seconds. The temperature at 1500 seconds can be estimated to be 958 K, 
according to the slope of temperature with time. This temperature difference of 77 
K may result from using a different emissivity for the heat pipe surface, neglecting 
natural convection in the numerical calculation, and additional heat loss from the 
inside surface and end of heat pipe through insulation. These effects increase with 
temperature. The leading edge test set was coated with a high emissivity cer ami c 
paint to improve heat rejection by radiation, but the emissivity was not specified. 
For numerical calculations an emissivity of 0.8 is used. The heat loss due to con- 
vection is estimated by using the approximation of a horizontal heated plate facing 
upward. Estimation shows that the operating temperature drop of 35 K may be 
possible due to convection. Increasing emissivity from 0.8 to 0.9 would reduce the 
outside temperature about 30 K. 

Figure 6.30 shows the temperature distributions at different wall locations. 
Near the stagnation point, a large temperature difference exists due to intensive 
heat input while for the rest of the heat pipe the difference is only about 15 K. The 
figure shows the length of the evaporator, which is not specified initially. Out to a 
length of about 15 cm, the temperature at the surface is greater than that of the 
interface, so that this region may be assumed to be the evaporator, and the rest of 
the length is the condenser section. Figures 6.31 and 32 show detailed temperature 
distributions obtained by numerical calculation, along the length at the external 


Time = 1450 Sec. 
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Figure 6.30. Temperature distribution at the different locations in wall. 
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surface and interface between shell and wick structure, respectively, for different 
times. 

Figure 6.33 shows the variation of the average vapor temperature during start- 
up. Vapor temperature increases corresponding to increasing heat flux. The vapor 
temperature increases very slowly because most of the heat added in the evaporator 
is extracted at the interface of the condenser to heat up the adjacent cold zone until 
the entire vapor space is in the continuum flow regime. Then the vapor temperature 
rises to approach steady state. Figure 6.34 shows the vapor flow dynamics at a time 
of 1140 seconds. Even though the sonic limit is not encountered, the velocity of the 
vapor is about 180 m/s. The temperature drop in the vapor space is about 60 K 
and pressure is not recovered due to the effect of friction at the interface, so that 
thermal resistance in the vapor space should be considered. 

Even though numerical results do not exactly match experimental results, due 
to the difference in configuration, material properties and boundary conditions, it 
is clear that the model does approximately predict the correct startup time and 
temperature distribution. It is important to note that heat inputs for the model 
were based on aerodynamic heating computations, which were only approximated 
in the experiments. 

To show the effectiveness of the heat pipe for cooling a leading edge, numerical 
calculations were executed for the same example problem, except that an adiabatic 
boundary condition is applied to the liquid-vapor interface during the entire startup. 
The development of large axial temperature gradients is similar for the cooled and 
uncooled leading edges, due to the free molecular condition in the entire vapor space 
during the initial transient heating. However, beyond this stage, peak temperatures 
are reduced, and temperature gradients disappear when the heat pipe is fully active. 
Figure 6.35 shows temperature distributions for both cases at a time of 1450 seconds. 


Temperature in Vapor Space 
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Figure 6.33. Average vapor temperature variation during startup. 





Figure 6.34. Axial variation of vapor temperature, pressure, velocity, and depsity 
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It is clear that the heat pipe is very effective for reducing peak temperatures near 
the leading edge and eliminating axial temperature gradients. 

6.4 Simulation of reentry heating environment 

After the model components were checked against published data, a spacecraft 
reentry heating environment [30] was simulated to show capability of the present 
model. The same rectangular cross-section and material of each component were 
used. Heat input at the stagnation point is much greater than that for the previous 
test, as shown in Figure 6.36, so that the heating distribution on the heat pipe 
surface is also much greater as shown in Figure 6.37. Maximum heat flux at the 
stagnation point is 390 kW/m 2 . A triangular element is used to fit the curved 
leading edge shape. Figure 6.38 shows the two-dimensional grid system used to 
represent the leading edge. In order to reduce computational time a relatively 
small number of elements and nodes was used. Heat input near the stagnation 
point is very large so that small size elements are used near the stagnation point 
and large size elements are employed for the rest of the heat pipe. The same 
numerical procedures as previously described were employed except that an initial 
time step of 5 seconds was used. 

Figure 6.39 shows temperature distributions along the heat pipe from the initial 
condition to 700 seconds. General behavior is similar to previous results, but due 
to high heat fluxes the continuum flow region is established early in the vapor space 
at 230 seconds. Also, as shown in Figures 6.40 and 41, most of the working fluid 
which is initially in the frozen state is melted at about 350 seconds and the entire 
vapor space is occupied by continuum flow at 460 seconds. As expected, less time 
is consumed for the entire heat pipe to become active. At a time of 500 seconds, 
the heat pipe is nearly isothermal with position, but the temperature increases 
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Figure 6.37. Heat flux distribution on a leading edge model for 
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uniformly with time until a maximum temperature of 1040 K is achieved at a time 
of 700 seconds. Figure 6.42 shows the temperature histories of the heat pipe at 
the stagnation point, end of the heat pipe, and a short distance away from the 
stagnation point. After the continuum flow region is established, it takes about 280 
seconds for the entire vapor space to reach continuum flow. Also Figures 6.43 and 
44 show temperature distributions along the length at an interface between shell 
and wick and at the liquid-vapor interface respectively, from initial time to 700 
seconds. Temperature distributions at these sections axe shown on Figure 6.45. A 
maximum temperature difference of 90 K is observed owing to the large heat input 
at the stagnation point while small differences exist for the rest of the heat pipe. 

Variation of the average vapor temperature as shown in Figure 6.46 is similar 
to that for low heat input. When the entire vapor space reaches the continuum 
flow regime, vapor temperature increases relatively rapidly and then again slows 
down to approach steady state. Figure 6.47 shows behavior of the vapor flow at 
503 seconds. A temperature drop of 50 K is observed along the length and pressure 
is not recovered. Since velocity is still large and heat extracted at the beginning of 
the condenser is small, the effect of friction at the interface is greater than that of 
inertia. Thus, a maximum velocity of the vapor does not occur at the exit of the 
evaporator but in the condenser. 

As expected, behavior of the heat pipe is similar to the previous results pre- 
sented except that a higher steady state temperature is achieved and startup is 
faster due to greater heat input. Figure 6.48 shows temperature distributions for a 
leading edge with heat pipe cooling and without cooling at 700 seconds into reentry. 
The heat pipe reduces peak temperature by the 600 K near the stagnation point 
and the leading edge becomes nearly isothermal. Again, it is obvious that the heat 
pipe is an effective device for cooling leading edges during reentry. 
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Figure 6.42. Temperature history of three locations in the heat pipe for reentry simulation. 
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Figure 6.43. Temperature distribution at interface between shell and capillary 
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Figure 6 . 44 . Temperature 


Time = 700 Sec. 
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Figure 6.46. Average vapor temperature variation during startup for reentry simulation 
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Figure 6.47. Axial variation of vapor temperature, pressure, velocity, and density 
at t = 503 seconds for reentry simulation. 
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Figure 6.48. Computed temperature distribution for uncooled and cooled leading edge 
of reentry simulation; t = 700 seconds. 
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CHAPTER VII 

CONCLUSIONS AND RECOMMENDATIONS 

7.1 Conclusions 

This study has been concerned with startup behavior of heat pipes from the 
frozen state. A mathematical model and associated finite element computational 
code were developed. Results were compared with published data and a reentry 
heating environment was simulated to show capability of the model developed. The 
model approximately predicts the correct startup time and temperature distribu- 
tion for the reentry problem. Numerical results computed for heat pipe cooled and 
uncooled leading edges show development of the large chordwise temperature gra- 
dients during the initial heating for both cases. After a part of the heat pipe is 
active, the effectiveness of using the heat pipe is shown clearly. In addition, the 
following conclusions are derived from numerical results. 

1. Temperature near the stagnation point increases rapidly and a large tem- 
perature gradient is observed due to the extremely small vapor density 
during the beginning of startup. Thus, critical design consideration should 
be given for this period. 

2. During the second phase of startup, the sonic vapor limit is encountered 
due to the large temperature gradient in the vapor space, and temperature 
increases slowly. 

3. When the sonic limit is not present, temperature first increases rapidly and 
then more slowly as steady state is approached. 

4. A large temperature difference exists at all times in the wall and capillary 
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structure near the stagnation point due to large heat input at this point. 

5. Heat pipe cooling of hypersonic wing structures greatly reduces tempera- 
ture gradients on the skin. 

6. Startup behavior of heat pipe cooled leading edges is similar for relatively 
large and relatively small heat inputs but large inputs cause startup to 
occur more quickly. 

7.2 Recommendations 

It is recommended that models and numerical techniques be improved so that 
more accurate predictions can be made and so that computer costs can be reduced. 

1. During startup, the working fluid is only partially melted so that some of 
the vapor may condense on the surface of the frozen working substance. 
Drying and rewetting might be incorporated into the model instead of 
assuming a saturated wick structure. 

2. For a better prediction of the vapor flow dynamics, one-dimensional, tran- 
sient and compressible equations might be developed to take into account 
supersonic flow of vapor. 

3. A better numerical scheme is needed to couple the effect of vapor flow 
dynamics to governing equations for the heat pipe shell and wick structure. 

4. To handle realistic operating conditions, the effects of gravity might be 
included in the mathematical model. 

5. The model should be made fully three dimensional. 


6. The computer program should be made user friendly. 


APPENDICES 



134 


APPENDIX A 

ELEMENT MATRICES 

A.l Interpolation function 

The spatial domain is discretized into three-node, linear, triangular elements 
as shown in Figure A.l, since an assemblage of triangles can always represent a 
two-dimensional domain of any shape. The shape functions are derived by using 
natural coordinates which range between zero and unity within the element, and 
whose variation between nodes is linear. Use of these coordinates is advantageous 
in evaluating the integrals in the element equations [55]. 

If Li, L 2 , and L 3 are selected as the natural coordinates, the location of the 
point P 0 within the element may be expressed by the following equations 


X = L\X 1 -|- L 2 X 2 -(- L 3 X 3 


(A.l) 


Y = L{Y 3 + L 2 Y 2 + L 3 Y 3 


(A.2) 


1 = L\ + L 2 + L 3 


(A3) 


From the equations above, natural coordinates are written in terms of cartesian 
coordinates. 


Li 


1 

2A< e > 


(dj —6] A + cjl ) 


(A4) 
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U = ^)(“2 + W + c 2 Y) (A.5) 


L> = ^j(Oi + b>X + c s Y) (A6) 


where 


ai=X 2 Y s -Y i Y 3 , b 1 = y 2 -y 3 , Cl =x 3 -x 2 


o 2 = x 3 yi-y 1 r 3 , &2 = n-^, C2 = x 1 -x 3 


a 3 = x a y 2 -y a y 2 , fe 3 = Fi-r 2 , c 3 = x 2 -x 3 

i JTi y 
2A< e > = i x 2 y 2 
i x 3 y 3 

These natural coordinates Li, L 2 , and i 3 are the linear interpolation functions 
for a triangle, that is, N{ = Li for the linear triangle. The interpolation functions 
are expressed as follows: 

Ni =2^) ( a i+ b iX + CiY) i = 1,2,3 (-A.7) 


(A. 8) 
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dN t b t 
dX ~ 2 A< e ) 
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dNi Ci 
dY ~ 2 A( e ) 


(A. 9) 


A.2 Two-dimensional element matrices 


Thermal conductivity and specific heat are assumed to be constant within an 
element and axe evaluated by using the average temperature. The [K c ] matrices are 
then expressed by 


where 


[K e 


■l. 


Ci 


— k( — 4- — Cj ^ 

\2A( e ) 2A< e ) 2A( e > 2A< C > ) 

= 4 (aw)’ (Wi + CiCi> 


(A.10) 




b\ + c\ 
b 2 bi + Cl c 2 
b 3 b\ + C 3 C 1 


b\ b 2 + ci C 2 
b\ + c\ 
b 3 b 2 + C 3 C 2 


£>1^3 + Ci C 3 

£>2 £*3 4- C 2 C 3 


£>7 4 - c 


For the element which is not involved in phase change, the [C] matrices are 
given by 


[C} = [ CNiNjdXdY 

JR ( •> 

= cf 

Jr(‘ 


Nf NjN 2 
N 2 N 1 Ni 
N 3 N 1 N 2 N 3 


N 1 N»\ 
iV 2 A f 3 
A 7 | / 


/ 1/6 1/12 1/12 \ 
= CA (e) 1/12 1/6 1/12 

\ 1/12 1/12 1/6 / 


(a.ll) 
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For the matrix contributions of an element having sides coincident with the 
external boundary of the solution domain, it is assumed that the convective heat 
transfer coefficient, h cr , heat flux, Q , and radiation coefficient, fi r , are constant 
along the side of am element. The shape functions for the boundary segment have 
the same expression of Equation(A.7). For the side of the element shown in Figure 
A.l between nodes 1 and 2, coincident with external boundary of the solution 
domain, the matrix contributions are written as follows: 



[K r }= f PrNiNjdS 

Ja* 


12 



(4.12) 


(4.13) 


where 1 n is the length of the side between nodes 1 nad 2 of the element. 
Also, the column vectors contributions from surface integrals are 



(4.14) 
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LISTING OF THE COMPUTER PROGRAM 


o o o 
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PROGRAM HPMAIN (INPUT , OUTPUT , TAPE3-INPUT , TAPEA-OUTPUT , HPDAT , 
5.TAPE5-HPDAT , JANG , TAPE6- JANG , DATA, TAPE7-DATA, RESTA, TAPE8-RESTA, 
S.DATAV , TAPE2-DATAV) 


C 
C 

c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c * 
c * 
c * 
c * 
c * 
c * 
c * 
c * 
c * 

c * 
* 


***************** * *************** ** * * * * * * * * * * * * * ******** * ******** * * * * * * 

* This program solves the startup and transient performance of heat * 
pipe with metallic working fluids by using a finite element method. 

The temperature is predicted from two-dimensional and transient 
heat conduction equation, which incorporates the effects of the 
phase change process in an expression for the volumetric heat cap- 
pacity by using the enthalpy method. The Galerkin weighted residual 
method is used to drive finite element formulations. 

The flow dynamics of the vapor are described by one-D, compress- 
ible and laminar momentum and energy equations. In one-D model, 
the variation of velocity at cross section, friction at the iiquid- 

* vapor interface, and quality of vapor are considered. 

Five one-D governing differential equations for vapor are solved 

by subroutine DVERK which uses Runge-Kutta method. The specified 
temperature, heat flux, convective, and radiation boundary 
conditions are applicable. This program can be used to solve pure 
conduction or phase change problem alone. 

Implicit or explicit time stepping schemes are used. Since 
explicit scheme is not self-s tatring method, implicit scheme is used* 
for first few time steps. For this purpose, set proper number for * 
variable NTS. 

Grid system can be generated by program HPGRNW. Input data file :: 
for program HPMAIN consists of output data file of HPGRNW and * 

general data which specifies some general conditions. 


To use IMSLIB Library : 

AT, IMSLIB/UN*LIBRARY 
LIBRARY, IMSLIB 

To obtain, compile and run program HPMAIN : 

GET, HPMAIN, HPDAT, RESTA (Get this RESTA only for restart) 
RWF 

FTN5,I*HPMAIN,L-0 
LGO, , 

(PROGRAM WILL PROMPT FOR INPUT DATA) 

To view results in files JANG, DATA and DATAV : 

RWF 

C, JANG (Temperature distribution for every time step) 
C,DATA(Temp. distribution for particalur time to plot) 
C, DATAV (Vapor temperature for every time step) 


* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 


C * 


Created ; Dec. 9, 1985 by Jong Hoon Jang 
Last update ; Dec. 05, 1987 by Jong Hoon Jang 



14 


142 


FORMAT ('HOW OFTEN DO YOU WANT TO PRINT, EX - 1,3..?') 
READ (3,*) NFREQ 


C 


C 

C 

C 

C 

C 

C 

c 

c 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


READ GENERAL DATA 

READ (5 , *) TEMP I , TMEL , TSTAR , DETP , HRSL , DIST .WIDTH 

TEMPI - INITIAL TEMPERATURE (TEMP.) [K] 

TMEL - MELTING TEMP [K] . 

TSTAR- TRANSIENT TEMP. OF VAPOR [K] 

DETP - TEMP. DIFFERENCE FROM TMELW 

HRSL - LATENT HEAT OF MELTING (ALWAYS POSITIVE) [J/KG] 

DIST - HEIGHT OF VAPOR SPACE [M] 

WIDTH - ELEMENT THICKNESS [M] 

TLENG - TOTAL LENGTH OF HEAT PIPE[M] 

READ (5 , *) BETA , DELT 1 , DELTP , MM , NLEM , NTS , NFLUX , NRAD 

BETA - i FOR FULLY IMPLICIT, .5 FOR CRANK-NIHOLSON 
DELT1 - TIME DIFFERENCE FOR FIRST TIME STEP[SEC.] 

DELT - TIME DIFFERENCE BETWEEN EACH TIME STEP [SEC.] 

MM - TOTAL NUMBER OF ITERATION FOR EACH TIME STEP 

MN - TOTAL NUMBER OF TIME STEP 

NLEM - 1 FOR LEMMON METHOD OF ENTHALPY 

2 FOR DEL GUIDICE METHOD OF ENTHALPY 
NFLUX - 1. HEAT FLUX DEPEND ON TIME AND CALL FLUXD 
NRAD - I RADIATION BOUNDARY CONDITION APPLIED. 

NTS - NUMBER OF STEPS OF IMPLICIT 

READ ELEMENT DATA PROVIDED FROM GRID2D 


C 

C 

C 

C 

C 

C 

C 

C 

C 

c 

c 

c 

c 


READ (5 , *) IOPT, LBW , NEL , NP , ICORD , IPROP , INRG.NDIM, IFIN 

IOPT - 0 FOR STEADY STATE, 1 FOR TRANSIENT 

LBW - BANDWIDTH QUANTITY 

NEL - NUMBER OF ELEMENTS 

NP - NUMBER OF NODES 

ICORD - 0 FOR RECTANGULAR COORDINATES 

1 FOR CYLINDRICAL COORDINATES 
IPROP * 0 FOR UNIFORM PROPERTIES THROUGHOUT 

1 FOR NON-UNIFORM PROPERTIES 
INRG - NUMBER OF REGION 

NDIM - 1 FOR ONE DIMENSIONAL PROBLEM 

2 FOR TWO DIMENSIONAL PROBLEM 

IFIN - 0, NOT A FIN PROBLEM 

1, A FIN PROBLEM 

READ (5 , *) (IRGN(I) ,NROW(l) ,NC0L(I) , 1-1, INRG) 


C 

C 


NCOL - NUMBER OF COLUMNS 
NROW - NUMBER OF ROW 
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I - 1 

INRG - INRG/ 2 
NCOLT = 0 
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C DECLARATION 

COMMON/MATG/ AMC(140,151) ,RM(140) ,THETA(140) ,FLD(140) , 

& AM (140, 140) , BM(140) ,CIJ(3,3) ,KIJ(3,3) 

4 /TEPS/ TEMPR, TEMPI, TMEL, DETP , DELTD, TEKPG1(140) ,THETAG(140) , 

4 THETAF (140) ,THETAO(140) ,THETTO(140) , THETAB (70) , TEMPD 

4 /PROC/ HRSL , ALPAR , DIST , B I , KREF , ROUR , SPECR , EMIS , HTC (140,3) , 

4 BOLT , VOLSPR , WIDTH , QE ( 1 40 , 3) .NRESIS ,RESIS , IX 

4 /ELMT/ NEL,NP,GND(140, 3) ,B (140, 3) ,C (140, 3) , EAREA(140) ,NPP1 , 

4 IES(140,2) ,ESL(140,3) , BETA, XC (140,3) , YC(140,3) ,DEF0, 

4 NE (140) , TLENG , PNOD (5,30) ,EBC(3,3) ,CM(3,3) ,GNDB(70) 

4 /INOUT/TEMPS (70) , TEMPV (70) , TEMPVS (70) , QEI (70) , VX (70) , 

4 GNDI (70) , VYl (70) ,VY2(70) ,VY3(70) ,ELSI(30) , 

4 VY4 (70) , VY5 (70) ,VY6(70) ,HFG,TEMPF (70) 

4 /SPLIN/TIME (3) .TIMER (20) ,QSTAG(20) ,XP1(25) , CPI (4,25) ,ND1, 

4 XP2 (60) ,CP2 (4,60) ,ND2,XP3(25) ,CP3(4,25) ,ND3,TS, 

4 XP4 (25) , CP4 (4,25) ,ND4,XP5(15) ,CP5(4,15) ,ND5,XP6(15) , 

4 CP6 (4, 15) ,XP7(15) ,CP7(4,15) ,ND7,XP8 (15) ,CP8 (4, 15) , 

4 ND6,ND8,XP9(15) ,CP9(4,15) ,ND9,XP0(15) ,C?0(4,15) ,ND0, 

4 ELSS (30) , TESL (30) 

4 /OUTPT/MO , IDAY , I YEAR .TITLE , CASE , NUMBER , NO , TEMPCC , HTCC , IOUT , 

4 T IMEN , DELT , DELT 1 , DELTP 

DIMENSION X (3) ,Y(3) ,RS(70) ,DUMG(140) ,DUMC(140) ,DUMQE(140,3) , 

4 IRGN (15) ,NCOL(15) ,NROW(15) ,DUMK(140) , TEST (140) , 

4 TEMP (3) , TEMPG (3) ,TEMPB(70) ,THETAX(140) ,TA(3,96) , 

4 AREA (30) , QT (70) ,QEIN(70) ,DTEMPS(70) , TEMPC (140 , 3) 

INTEGER TS , TSM1 , TSM2 , P , GND , GNDB , GNDI , IOUT , FIRTS , FMl , FP1 , PNOD 
REAL KI J , KLP , KSN , KFN , KLN , KREF , HRSL , HTC , KEE , MACH , INTVL, MASW , MASWN , 
4 JCOB , JCOBI 

C READ IN DATA TO CREATE COVER PAGE 


WRITE(4,5) 

5 FORMAT (' INPUT CASE NUMBER?') 

READ (3 , *) NO 

WRITE (4,6) 

6 FORMAT (' INPUT DATE AS 12 20 1987') 

7 FORMAT (13 , 13 , 15) . 

READ (3, 7) MO, IDAY, I YEAR 
WRITE(4,8) 

8 FORMAT (’ INPUT TITLE WITHIN 27 CHARACTERS') 

9 FORMAT (3A9) 

READ (3, 9) TITLE, CASE, NUMBER 

10 FORMAT (Al) 

WRITE (4, 11) 

11 FORMAT ("RESTART FROM LAST RUN ? (Y/N) ") 

READ (3 , 10) YESNO 

WRITE(4, 12) 

12 FORMAT ( ’ TIME INTERVAL (SEC .)?' ) 

READ (3,*) DELT 

WRITE(4, 13) 

13 FORMAT (’HOW MANY TIME STEP DO YOU WANT TO RUN?') 
READ (3 , *) MN 

WRITE (4,14) . 


ooooo oononon oooo 
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NROWT - NR0W(1) + NR0W(INRG+1) - 1 
30 NCOLT - NCOLT + NCOL(l) 

IF (I .LT. I NRG) THEN 
1*1 + 1 
GO TO 30 
END IF 

NCOLT - NCOLT - INRG + 1 
DO 100 I - l.NEL 

READ (5 , *) NE (I) , (GND(I.L) ,L-1,3) , (XC(I,L) ,YC(I,L) ,L-1,3) 

NE “ NUMBER OF ELEMENT 
GND - GLOBLE NUMBER OF NODES 
XC - X-COORDINATE OF NODES [M] 

YC - Y-COORDINATE OF NODES [M] 

READ (5,*) DUMK(I) ,DUMG(I) ,DUMC(l) , (IES(I.L) , 

& HTC (I , IES (I , L) ) , TEMPC (I , IES (I, L) ) , L*1 , 2) , 

& (IES(I.L) , QE (I , IES (I , L) ) , L*1 , 2) 

DUMK - THERMAL CONDUCTIVITY [W/ (M*K) ] 

DUMG - HEAT GENERATION PER UNIT VOLUME 
DUMC - SPECIFIC HEAT/DENSITY PRODUCT 
HTC - HEAT TRANSFER COEFFICIENT [W/ (M**2*K) 3 
IES - NUMBER OF SIDE IN EACH ELEMENT 
QE - HEAT FLUX ON BOUNDARY SURFACES [W/M**2] 

TEMPC - REFERENCE TEMP. FOR CONVECTIVE HEAT [K] 

100 CONTINUE 
1ST - 1 
IFN - 6 

40 READ (5 , *) (GNDB(I) ,1 - 1ST, IFN) , (TEMPB(I) ,I-IST,IFN) 

GNDI - GLOBLE NUMBER OF NODE WHERE IS LOCATED ON 
THE LIQUID-VAPOR INTERFACE. 

GNDB - GLOBLE NODE NUMBER OF NODE WHERE IS ON BOUNDARY 
JFNT - TOTAL NUMBER OF GNDI 
TEMPB - SPECIFIED BOUNDARY TEMP. 

IF (GNDB (IFN) .GT. 0) THEN 
1ST - 1ST + 6 
IFN - IFN + 6 
GO TO 40 
END IF 

C READ GLOBAL NODE NUMBER ON LIQUID- VAPOR INTERFACE 

JST - 1 
JFN - 6 

■50 READ (5 , *) (GNDI (I), I - JST, JFN) 

IF (GNDI (JFN) . GT. 0) THEN 
JST - JST + 6 
JFN = JFN + 6 
GO TO 50 
END IF 
JFNT - 0 
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60 


DO 60 I - l.JFN 

IF (GNDI (I) .GT. 0) THEN 
JFNT - JFNT + 1 
END IF 
CONTINUE 


C READ DATA FOR RADIATION BOUNDARY CONTITIONS 

IF (NRAD .EQ. 1) THEN 
READ (5 , *) TEMPR.EMIS 
END IF 


C TEMPR - REFERENCE TEMP. FOR RADIATION [K] 

C EMIS - EMISSIVITY 

C READ SCALE FACTOR AND HEAT FLUX AT STAGNATION 


15 FORMAT (A20) 

IF (NFLUX .EQ. 1) THEN 
READ (5, 15) INPUTl 

READ (5 , *) ND1, (XP1(I) , CPI (1 , I) , 1-1 ,ND1) 

READ (5 , *) CPI (2, 1) ,CP1 (2, ND1) 

CALL SPLINE (ND1-1.XP1, CPI) 

CALL CALCF(ND1-1,XP1,CP1) 

READ (5, 15) INPUT2 

READ (5 , *) ND2, (XP2 (I) ,CP2(1, 1) , I-1.ND2) 

READ (5 , *) CP2 (2, I) ,CP2 (2,ND2) 

CALL SPLINE (ND2-1,XP2,CP2) 

CALL CALCF (ND2- 1 , XP 2 , CP 2 ) 

END IF 

C READ CONDUCTIVITY AND SPECIFIC HEAT OF STAINLESS STEEL 

READ (5, 15) INPUT3 

READ (5 , *) ND3 , (XP3 (I) , CP3 ( 1 , I) , 1-1 , ND3) 

READ (5,*) CP3 (2,1) , CP3 (2 , ND3) 

CALL SPLINE (ND3-I,XP3,CP3) 

CALL CALCF (ND3-1.XP3.CP3) 

READ (5, 15) INPUT4 

READ (5,*) ND4, (XP4(l) ,CP4(1 , 1) , 1-1 ,ND4) 

READ (5,*) CP4(2, 1) ,CP4(2,ND4) 

CALL SPLINE (ND4-1.XP4.CP4) 

CALL CALCF (ND4-1 ,XP4,CP4) 


C NDl.2,3 - NUMBER OF DATA POINTS 

C XP 1 , 2 , 3 - LOCATION IN X - AXIS 

C CP 1,2, 3 - COEFFICIENT OF INTERPOLATION EQUATION 

C CONDUCTIVITY .DENSITY AND SPECIFIC HEAT OF SOLID SODIUM 


READ (5, 15) INPUT5 

READ (5,*) ND5, (XP5 (i) ,CP5(1,I) .I-1.ND5) 
READ (5,*) CP5 (2,1), CP5 (2 , ND5) 

CALL SPLINE (ND5-1.XP5.CP5) 
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CALL CALCF(ND5-1,XP5,CP5) 

READ (5, 15) INPUT6 

READ (5,*) ND6 , (XP6(I) ,CP6(1,I) ,I-1,ND6) 

READ (5 , *) CP6 (2, 1) ,CP6 (2,ND6) 

CALL SPLINE (ND6-1,XP6,CP6) 

CALL CALCF(ND6-1,XP6,CP6) 

READ (5, 15) INPUT7 

READ (5 , *) ND7 , (XP7 (I) , CP7 (1 , 1) , 1-1 , ND7) 

READ (5,*) CP7 (2 , 1) , CP7 (2 , ND7) 

CALL SPLINE (ND7-1.XP7.CP7) 

CALL CALCF (ND7- 1 , XP 7 , CP 7) 

....CONDUCTIVITY .DENSITY AND SPECIFIC HEAT OF LIQUID SODIUM 
READ (5, 15) INPUT8 

READ (5,*) ND8, (XP8(I) ,CP8(1,I) .I-1.ND8) 

READ (5,*) CP8 (2, 1) ,CP8 (2 ,ND8) 

CALL SPLINE (ND8-1.XP8.CP8) 

CALL CALCF (ND8-1.XP8.CP8) 

READ (5 , 15) INPUT9 

READ (5.*) ND9 , (XP9(I) ,CP9 (1 , I) , I-1.ND9) 

READ (5,*) CP9 (2,1) , CP9 (2 , ND9) 

CALL SPLINE (ND9-1.XP9.CP9) 

CALL CALCF (ND9-1.XP9.CP9) 

READ (5, 15) INPUTO 

READ (5,*) NDO, (XPO(I) ,C?0(1 , I) . I-l.NDO) 

READ (5,*) CP0(2, 1) ,CP0(2,ND0) 

CALL SPLINE(NDO-l.XPO.CPO) 

CALL CALCF (NDO- l.XPO, CP 0) 


C READ NODAL POINT NUMBER TO USE PRINTING TEMP 

READ (5,*) ((PNOD(I, J) , J-l.NCOLT) ,1-l.NROWT) 

C PRINT OUT DATA FOR PLOT 

IF (YESNO .EQ. "N") THEN 

WRITE (7.*) ((GND(I.L) , L-1,3) , I-l.NEL) 

WRITE (7,*) (CXCCl, L).YC(I.L), L-1,3) ,1-1, NEL) 
END IF 


C SET FIRST GUESS 


DO 110 I- 1ST , IFN 
IF (GNDB (I) .GT. -1) THEN 
TEMPGl(GNDBd)) - TEMPB(I) 
END IF 

110 CONTINUE 

DO 120 I - 1,NP 
DO 125 J - 1ST, IFN 

IF (I .NE. GNDB (J) ) THEN 
TEMPG1 (I) = TEMPI 
END IF 

125 CONTINUE 
120 CONTINUE 
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C SOME CONSTANT PARAMETER 

JNT - NUMBER OF NODES ON INTERFACE FOR PRESENT TIME STEP. 

TS - NUMBER OF TIME STEPS 


BOLT « 5.67E-8 
FM1 - 1 
FIRTS - 1 
FP1 - 1 
KLP - 72.3 
KREF - 141. 

VOLSPR - 1.109E6 
ALPAR - KREF/ (VOLSPR) 

DEFO - ALPAR*DELT1/DIST**2 
IOUT - 1 
JK - 1 
JNT - 1 
IX - 1 

NPPl - NP + 1 
NRESIS - 1 
TIMEN - DELT1 
■ TIME (I) - DELT1 
TS - 0 

TEMPV (2) - 0. 

C CALCULATE BASIC ELEMENT GEOMETRIC MATRIX 

DO 200 M - l.NEL 
DO 220 I - 1,3 
X(I) - 0. 

Y (I) - 0. 

X(I) - X(I) + XC(M,I)/DIST 
Y (I) - Y(I) + YC(M,I)/DIST 
220 CONTINUE 

CALL ELMDM (X , Y , M) 

200 CONTINUE 

C CALCULATE LENGTH OF SIDE OF ELEMENT AT INTERFACE 

K = 2 
KM1 - 1 

226 DO 225 M » 1 ,NEL 

IF (GND (M, 1) .EQ. GNDI (KM1) .AND. GND (M, 2) .EQ. GNDI (K) ) THEN 
ELSI(KMl) - SQRT ( (XC (M, 1) - XC ( M,2))**2 + 

& (YC (M, 1) - YC (M, 2) ) **2) 

K - K + 1 
KM1 = K - 1 

IF (K .LE. JFNT-1) THEN 
GO TO 226 
END IF 
END IF 
225 CONTINUE 
K = 2 
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KM1 “ 1 
TLENG - 0.0 

228 DO 227 8 - l.NEL 

IF(GND(M,3) .EQ. PN0D(1,KM1) .AND. GND(M,2) .EQ. 
& PN0D(1,K)) THEN 

ELSS(KMl) - SQRT ( (XC (M, 2) - XC(M,3))**2 + 

& (YC(M,2) - YC (M, 3) ) **2) 

TLENG - TLENG + ELSS (KMl) 

IF (KM1 .EQ. 1) THEN 
TESL(l) - ELSS(l) 

ELSE 

TESL(KMl) - TESLCKMl - 1) + ELSS (KMl) 

END IF 
K - K + 1 
KM1 - K - 1 

IF (K .LE. JFNT-1) THEN 
GO TO 228 
END IF 
END IF 
227 CONTINUE 

WRITE(7,*) TLENG, (ELSS (I) ,1-1, JFNT-1) 

WRITE (7,*) (TESL(I) ,1-1, JFNT-1) 

WRITE (7,*) (ELSI (I), 1-1, JFNT-1) 

IF(YESNO .EQ. "Y") THEN 
GO TO 545 
END IF 


q ********************* * * * * * * * * ******************** 

C TO EVALUATE TEMPERATURE AT FIRST FEW TIME STEPS 
C BY USING IMPLICIT METHOD (BETA = 1...5) 
c ************************************************* 

C SET DIMENSIONLESS BOUNDARY TEMP 

DO 230 I - l.IFN 

THETAB (I) - (TEMPB(I) - TMEL) / (TMEL - TEMPI) 

230 CONTINUE 

C SET DIMENSIONLESS GUESS TEMP 

DO 250 I - 1 , NP 

THETAG(I) - (TEMPG1 (I) - TMEL) /(TMEL - TEMPI) 

250- CONTINUE 
450 TS - TS + 1 

IF (TS .EQ. 2) THEN 

TIME (2) - TIME(l) + DELT 
ELSE IF (TS .EQ. 3) THEN 
TIME (3) - TIME (2) + DELT 
ELSE IF (TS .GT. 3) THEN 
IIME(I) = TIME (2) 

TIME (2) - TIME (3) 

TIME (3) = TIME (3) + DELT 
END IF 

DO 260 I - 1,NP 

TEMPG1 (I) = THETAG (I) * (TMEL - TEMPI) * TMEL 

ORIGINAL PAGE IS 
OF POOR QUALITY 


ORIGINAL PAGE IS 
OF POOR QUALITY 


260 CONTINUE 
430 P - 0 

400 P - P + 1 

C SET ARRAYS ZERO 

DO 310 I - 1 , NP 
THETA (I) - 0. 

310 CONTINUE 

DO 320 I - 1 , NP 
FLD(I) - 0. 

BM(I) - 0. 

RM(I) - 0. 

DO 325 J - 1 ,NP 
AM(l,J) - 0. 

325 CONTINUE 
320 CONTINUE 

DO 330 I - 1,NP 
DO 335 J - 1,NPP1 
AMC(I, J) - 0. 

335 CONTINUE 

330 CONTINUE 

DO 340 M - l.NEL 

DO 345 I - 1,3 

TEMPG(l) - THETAC (GND (M, I) ) * (TMEL - TEMPI) + TMEL 

C IMPLICIT TIME STEPPING SCHEME. 

IF (BETA .EQ. 1.) THEN 
TEMP (I) - TEMPG(I) 

ELSE 

TEMP (I) - (TEMPG(I) + TEMPG1 (I) ) *BETA 
END IF 
345 CONTINUE 

C TO ASSUMBLE ELEMENT MATRIX INTO GLOBAL SYSTEM MATRIX [AMC] 

CALL GLOBMAX (TEMP , TEMPC, M , DUMK , NLEX , NFLUX , NR AD , NTS , P , I FN) 
340 CONTINUE 


C APPLY B. C.'S OF FIRST KIND IF APPLICABLE 

CALL FIRSTBC (IFN) 

C TO EVALUATE RESIDUAL AND THETAG(I) FOR NEXT ITERATION 

DO 350 I - 1,NP 

THETA (I) - THETAG (I) 

THETAG(I) - 0. 


350 CONTINUE 

CALL RESIDUL(RMAX) 

CALL CHLSKY (TS , NTS , THETAX) 

DO 355 I = 1 . NP 

THETAG (I) = THETAX (I) + THETA (I) 
355 CONTINUE 



c 

c 

c 

c 

c 

c 

c 


*********************************************************** ******* 


* CHECK CONVERGENCE CRITERIOR ON DIFFERENCE BETWEEN NEW THETA AND * 
’•OLD THETA AT EVERY NODES. IF ONE OF DIFFERENCES IS GREATER THAN * 
’•TOLERANCE, THEN ITERATED BY NEWTON-RAPHSON METHOD. IF CRITERION * 
*IS MET, THETA (I) ARE SUBSTITUTED INTO EQUATION TO OBTAIN * 

*THETA(I) AT THE NEXT TIME STEP. * 

* **** **** ********* ***************** ********** *** * ******* ****** * *** 


IF(RMAX .LT. .001 .OR. P .GT. MM) THEN 
DO 370 I - 1,NP 
IF (TS .EQ. 1) THEN 

TA(1,I) - THETAG (I) * (TMEL - TEMPI) + TMEL 
ELSE IF (TS .EQ. 2) THEN 

TA(2,I) - THETAG (I)* (TMEL - TEMPI) + TMEL 
ELSE IF (TS .EQ. 3) THEN 

TA(3,I) - THETAG (I)* (TMEL - TEMPI) + TMEL 
ELSE 

TA(I.I) - TA(2,I) 

TA(2, 1) - TA(3,I) 

TA(3, I) - THETAG (I)* (TMEL - TEMPI) + TMEL 
END IF • 

370 CONTINUE 

NPRINT - JK * NFREQ 
IF (TS .EQ. NPRINT) THEN 

CALL MAINOUT (TA , TSTAR , NTS , MN , NCOLT , NROWT) 

JK - JK + 1 
END IF 

DEFO - (ALPAR/DIST**2) *DELT 
IF (TS .LT. NTS) THEN 
GO TO 450 
ELSE 

IF (NTS .EQ. MN) THEN 
GO TO 540 
ELSE 

GO TO 550 
END IF 
END IF 
ELSE 

GO TO 400 
END IF 

540 CONTINUE 
GO TO 2 


Q * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ****** * * * * * * * * * * * * * * * * * 

C TO EVALUATE TEMPERATURE BY USING EXPLICIT METHOD 
C (DUPONT) WITH TEMPERATURES FOR TWO PREVIOUS TIME STEPS 

Q ****** ************************************************** 


545 IF(YESNO .EQ. "Y") THEN 

CALL DATAIN (TA, JK , JNT , TIMEN , FIRTS , TVU , QT , TEMPD) 
FMl - FIRTS - 1 
FP1 = FIRTS 
END IF 
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C SET ARRAYS ZERO 

550 TS - TS + 1 

IF (TS .EQ. 2) THEN 

TIME (2) - TIME (1) + DELT 
ELSE IF (TS .EQ. 3) THEN 
TIME (3) - TIME (2) + DELT 
ELSE IF (TS .GT. 3) THEN 
TIME(l) - TIME (2) 

TIME (2) - TIME (3) 

TIME (3) - TIME (3) + DELT 
END IF 

530 DO 500 I - 1,NF 
FLD(I) - 0. 

BM(I) - 0. 

DO 510 J-l.NFPl 
AMC(I.J) - 0.0 
510 CONTINUE 
500 CONTINUE 

C SET INITIAL TEMPERATURE FROM TWO PREVIOUS TIME STEPS 

DO 520 I - 1 ,NP 
I*F (TS .GT. 3) THEN 

THETTO(I) - (TA(2,I) - TMEL) / (TMEL - TEMPI) 

THETAO (I) - (TA(3 , I) - TMEL) / (TMEL - TEMPI) 

ELSE 

THETTO (I) - (TA(1 , I) - TMEL) / (TMEL - TEMPI) 

THETAO (I) - (TA(2, I) - TMEL) / (TMEL - TEMPI) 

END IF 

DO 525 J - 1 ,NP 
AM(I,J) - 0. 

525 CONTINUE 

520 CONTINUE 

DO 600 M - l.NEL 

DO 610 I - 1,3 

TEMP (I) - THETAO (GND(M, I))* (TMEL - TEMPI) + TMEL 
610 CONTINUE 

C .... TO ASSUMBLE ELEMENT MATRIX INTO GLOBAL SYSTEM MATRIX [AMC] 

CALL GLOBMAX (TEMP , TEMPC , M.DUMK , NLEM , NFLUX , NRAD, NTS , P , IFN) 
600 CONTINUE 

C APPLY B.C.'S OF FIRST KIND IF APPLICABLE 

CALL FIRSTBC (IFN) 

C SOLVE MATRIX BY USING CHOLESKY 1 S DECOMPOSITION 

CALL CHLSKY (TS , NTS , THETAX) 

DO 646 I = 1 ,NP 

THETA (I) = THETAX (I) 

646 CONTINUE 
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C CHECK TEMPERATURES ON LIQUID-VAPOR INTERFACE 

C GNDI(l) - GLOBAL NODE NUMBER AT FIRST NODE ON INTERFACE 

C GNDI(JNT) - LAST GLOBAL NODE NUMBER WHOSE TEMPERATURE IS 

C GREATER THAN TSTAR 

C TEMPS - LIQUID TEMP. AT LIQUID-VAPOR INTERFACE 

C TEMPV - TEMP. IN VAPOR SPACE 

C TEMP VS - VAPOR TEMP. AT LIQUID-VAPOR INTERFACE 

DO 660 I - 1 ,NP 
IF (TS .GT. 3) THEN 
TA (1 , 1) - TA(2,I) 

TA(2, 1) - TA(3,I) 

TA(3, I) - THETA(l)*(TMEL - TEMPI) + TMEL 
ELSE 

TA(1,I) - TA(1 , 1) 

TA(2, 1) - TA(2, 1) 

TA (3 , I) - THETA (I)* (TMEL - TEMPI) + TMEL 
END IF 
660 CONTINUE 

NPRINT - JK * NFREQ 
IF (TS .EQ. NPRINT) THEN 

CALL MAINOUT (TA , TSTAR , NTS , MN , NCOLT , NROWT) 

JK - JK + 1 
END IF 

DO 665 I - l.JFNT 
TEMPS (I) - 0. 

665 CONTINUE 

DO 680 I - 1,2 
JJ - GNDI(I) 

TEMPS (I) - THETA (JJ) * (TMEL - TEMPI) + TMEL 
680 CONTINUE 

C ADIABATIC BOUNDARY CONDITION AT LIQUID-VAPOR INTERFACE 

IF (TEMPS (2) .LT. TSTAR) THEN 
DEFO - (ALPAR/DIST**2) *DELT 
IF (TS .GE. MN) THEN 
GO TO 2000 
ELSE 

DO 910 J - l.JFNT 
DO 911 I - l.NEL 

IF(GNDI(J) .EQ. GND(I , 1) .AND. GNDI(J+1) .EQ. GND(I,2)) 
& THEN 

IES(I.I) - 1 
QE(I,IES(I,1)) - o. 

END IF 

911 CONTINUE 

910 CONTINUE 

IF (TEMPS (JFNT) .GT. TSTAR) THEN 
GO TO 2000 
ELSE 

GO TO 550 
END IF 
END IF 
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ELSE IF (TEMPS (1) .GT. TSTAR .AND. TEMPS (2) .GT. TSTAR) THEN 

...NO LONGER ADIABATIC CONDITION AT LIQUID-VAPOR INTERFACE... 
.HOW MANY NODES ON INTERFACE HAVE GREATER TEMP. THAN TSTAR... 

DO 915 I - l.JFNT 
TEMPS (I) - 0. 

915 CONTINUE 

J * 1 

920 JG - GNDI(J) 

TEMPS (J) - THETA (JG) * (TMEL - TEMPI) + TMEL 
IF (TEMPS (J) .GT. TSTAR) THEN 
J - J + 1 

IF(J .LT.JFNT) THEN 
GO TO 920 
ELSE ' 

JNT - JFNT - 1 
END IF 
ELSE 

JNT - J 

TEMPSN - TEMPS (J) 

END IF 

C CALCULATE UNIFORM TEMP. FOR VAPOR SPACE 

TEMPS (JFNT) - THETA (GNDI (JFNT) ) * (TMEL-TEMPI) +TMEL 
DO 922 J - l.JFNT 

DO 923 I - l.NEL 

IF (GNDI (J) .EQ. GND(I.l) .AND. GNDI(J+1) .EQ. GND(I,2)) 
& THEN ‘ 

IES(I.l) - 1 

QE(I , IES (1 , 1) ) - 0. 

END IF 

923 CONTINUE 
922 CONTINUE 

IF (NRESIS .EQ. 2) THEN 

CALL COUPLE ( JNT, QT.TVU, JFNT) 

GO TO 2500 
END IF 

IF (FIRTS .EQ. 1) THEN 

TVU - (TEMPS (1) + TSTAR) /2. 

2100 RESD - 0. 

TELSI - 0. 

RESDC - 4.5E6*1.4*.021*2.29E11 
DO 924 I - 1, JNT 

RESDI - ELSI (I) *10. ** (-5567. /TEMPS (I) ) /TEMPS (I) 

RESDV - ELSI (I) *10. **(-5567. /TVU) /TVU 
RESD - RESD + RESDC* (RESDI - RESDV) 

TELSI - TELSI + ELSI (I) 

924 CONTINUE 

JCOB - (1 - 128 18 . 5/TVU) *TELSI 
JCOB - JCOB * 10. **(-5567. /TVU) 

JCOB = RESDC* JCOB/ (TVU**2) 

TVO - TVU 

TVU = TVO - RESD/ JCOB 
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IF(ABS(RESD) .LT. .001) THEN 

CALL INTFLUX(JFNT, JNT, TVU, TEMP V, TEMPS, QEI) 

ELSE 

GO TO 2100 
END IF 
ELSE 

2200 RESD - 0. 

TELSI - 0. 

RESDC - 4.5E6*1.4*.021*2.29E11 
DO 928 I - l.JNT 

RESDI - ELSI (I) *10. **(-5567. /TEMPS (I)) /TEMPS (I) 

RESDV - ELSI (I) *10.** (-5567. /TVU) /TVU 
RESD - RESD + RESDC* (RESDI - RESDV) 

TELSI - TELSI + ELSI(l) 

928 CONTINUE 

JCOB - (1 - 12818. 5/TVU)*TELSI 
JCOB - JCOB * 10. **(-5567. /TVU) 

JCOB - RESDC* JCOB/ (TVU**2) 

TVO - TVU 

TVU - TVO - RESD/ JCOB 
IF(ABS(RESD) .LT. .001) THEN 

C TO ESTIMATE THE SONIC LIMIT 

CALL INTFLUX (JFNT , JNT , TVU , TEMPV , TEMPS , QEI) 

VMAX = SQRT(1. 4*8314. *TVU/23.) 

PSATT = (2. 29E1 1/SQRT (TVU) ) *10. ** (-5567/TVU) 

DENS = 2.766E-3*PSATT/TVU 

HFGS - 182. *(25474. 93 - .9935*TVU) 

QMAX - VMAX*HFGS*DENS*DIST*WIDTH/2.2 
QTAL - 0. 

DO 930 I - l.JNT 
AREA (I) - 0. 

IF (QEI (I) .GT. 0.) THEN 

QTAL - QTAL + QEI (I) *ELSI (I) *WIDTH 
END IF 

930 CONTINUE 

QMAXT - QMAX*. 50 

.IF AMOUNT OF HEAT TRANSFER ON EVAPORATOR IS GREATER 

.THAN THE SONIC LIMIT, IT IS ASSUMED THAT THE SONIC 

.LIMIT IS ACTUAL HEAT TRANSFER. CALCULATE NEW VAPOR 

.TEMPERATURE BASE ON SONIC LIMIT 

IF (QMAX .LE. QTAL .AND. JNT .GE. 4) THEN 
2300 RESD - 0 

JCOB - 0. 

DO 931 I - 2, JNT 

IF (TEMPS (I) .GT. TVU) THEN 
JNTL = I 

ELSE IF (TEMPS (I) .LT. TVU .AND. TEMPS(l-l) .GT. TVU) 
4 THEN 

XII = (TEMPS (JNTL) - TVU) *ELS I (I- 1) 

XII = XII/ (TEMPS (JNTL) -TEMPS (I)) 

END IF 
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931 CONTINUE 

RSC1 - 2. 29E1 1*1 . 5*. 021*4 . 5E6 
RSC2 - 6 . 34E8*4 . 5E6* 10 . 26*DIST 
IC - 0 

DO 932 I - 1 , JNTL 

IF(QEI(I) .GT. 0.) THEN 
RSC12 - RSC1*ELSI (I) 

IF (I .NE. JNTL) THEN 

RESDI - 10.** (-5567/TEMPS (I)) /TEMPS (I) *ELSI (I) *RSC1 
RESDV - 10. ** (-5567/TVU) /TVU*RSC12 
JCOBI - RSC1*ELSI (I) 

ELSE 

RESDI - 10.** (-5567/TEMPS (I)) /TEMPS (I) *XII*RSC1 
RESDV - 10. ** (-5567/TVU) /TVU* (RSC1*XII+RSC2) 

JCOBI - RSC1*XII 
END IF 

RESD - RESD + RESDI - RESDV 
JCOB - JCOB + JCOBI 
ELSE 

IC « IC + 1 
END IF 

932 CONTINUE 

JLM1 - JNTL - IC 

JCOB -(JC0B+RSC2) *(1-12818. 5/TVU) 

JCOB - JCOB*10.** (-5567. /TVU) / (TVU**2) 

TVO - TVU 

TVU - TVO - RESD/ JCOB 
IF (ABS (RESD) .LT. .1) THEN 


C CALCULATE HEAT FLUX DISTRIBUTION ON CONDENSER SECTION 

AREAT - 0. 

DO 933 I - 2, JNT 

IF (TEMPS (I) .LT. TVU .AND. TEMPS (1-1) .GT. TVU) 

& THEN 


JNTM - I 

AREA (I) - TEMPS (I) - TVU 

AREA (I) - AREA (I) + (TEMPS (I) - TEMPS (1+1) ) /2. 
AREA(I) - AREA(l) *ELSI (I) 

ELSE IF (TEMPS (I) .LT. TVU) THEN 
AREA (I) - TEMPS (I) - TVU 

AREA(I) - AREA (I) + (TEMPS (I) - TEMPS (1+1) ) /2. 
AREA (I) - AREA (I) *ELSI (I) 

END IF 

AREAT - AREAT + AREA (I) 

933 CONTINUE 

JNTL - JNTM - 1 

CALL INTFLUX (JFNT, JNTL , TVU, TEMP V , TEMPS , QEI) 

VMAX - SQRT( 1.4*83 14. *TVU/23.) 

PSATT - (2. 29E1 1/SQRT (TVU) ) *10. ** (-5567/TVU) 

DENS =* 2.766E-3*PSATT/TVU 

HFGS = 182. *(25474. 93 - .9935*TVU) 

QMAX - VMAX*HFGS*DENS*DIST*WIDTH/2.2 
QTAL » 0. 

DO 935 I - 1 , JNTL 



n n 


156 


IF(QEI(I) .GT. 0.) THEN 
IF (I .LT. JNTL) THEN 

QTAL - QTAL + QEI (I) *ELSI (I) *WIDTH 
ELSE 

XII - (TEMPS (JNTL) - TVU)*ELSI(I) 

XII - XII/ (TEMPS (JNTL) -TEMPS (JNTM)) 

QTAL - QTAL + QEI (I) *XII*WIDTH 
END IF 
END IF 

935 CONTINUE 

DO 936 I - JNTM, JNT 
IF (I .EQ. JNTM) THEN 

QEI (I) - - AREA (I) * QMAX/AREAT/ELSI (i) /WIDTH 
IF(TS .GT. 7000) THEN 
END IF 
ELSE 

QEI (I) - - AREA(I) *QMAX/AREAT/ELSI (I) /WIDTH 
END IF 

936 CONTINUE 
NRESIS - 1 

ELSE 

GO TO 2300 
END IF 

ELSE IF (QTAL .LE. QMAXT .AND. JNT .GE. JFNT-1) THEN 


ARTIFICIAL THERMAL RESISTANCE, WHICH CALCULATED FROM 
VAPOR FLOW, IS IMPOSED AT LIQUID-VAPOR INTERFACE 


2500 


& 


I 


9A6 

945 


CALL COUP LE( JNT, QT.TVU, JFNT) 

NRESIS » 2 
END IF 
ELSE 

GO TO 2200 
END IF 
END IF 

IF (TS .EQ. NPRINT) THEN 
WRITE (6, 790) TVU 
WRITE (2, 775) TVU 
END IF 

DO 945 J - l.JNT 
DO 946 I - 1 ,NEL 

IF(GNDI(J) .EQ. GND(I,1) .AND. GNDI(J+1) .EQ. 
)) THEN 

IES(I,1) - 1 

QE(I,IES(I,1)) = QE(I,IES(I,1)) - QEI ( J) 

END IF 
CONTINUE 
CONTINUE 

FIRTS = FIRTS + 1 
DEFO = (ALPAR/DIST**2)*DELT 
IF (TS .GE. MN) THEN 
GO TO 2000 
ELSE 

GO TO 550 
END IF 


GND (1,2 
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ELSE 

DEFO - (ALPAR/DIST**2)*DELT 
IF (TS .GE. MN) THEN 
GO TO 2000 
ELSE 

GO TO 550 
END IF 
END IF 

C PRINT UNIFORM VAPOR TEMPERATURE 

2000 CONTINUE 
2 FIRTS1 - FIRTS - 1. 

775 FORMAT (F8. 2) 

790 FORMAT C/2X, 'TEMP V - ' ,F8.3, 2X, ' [K] ') 

C PRINT DATA FOR RESTARTING FROM LAST RUN 

CALL DATAOUT (TA , JK , JNT , TIMEN , FIRTS , TVU , QT , TEMPD) 

STOP 

END 


SUBROUTINE ELMDM(X, Y,M) 

q ****************************** ************************ 

C * THIS SUBROUTINE CALCULATES ELEMENT DATA SUCH * 

C * AS AREA, B (I) B ( J) + C(I)C(J), B(I)C(J) - B(J)C(l) * 

C * C1(I,J) AND CM(I, J) , AND FORMS MATRICES S622, S262 * 

C * AND S226 * 

q ****************************************************** 

COMMON/PROC/ HRSL , ALP AR , DI ST , B I , KREF , ROUR , SPECR , EMI S , HTC ( 1 40 , 3 ) , 

& BOLT, VOLSPR, WIDTH, QE(140, 3) ,NRESIS,RESIS, IX 

& /ELMT/ NEL,NP,GND(140,3) ,B(140,3) ,C(140,3) ,EAREA(140) ,NPP1, 

& IES (140, 2) ,ESL(140,3) ,BETA,XC (140 , 3) ,YC(140,3) , DEFO , 

& NE (140) , TLENG,PNOD(5 ,30) ,EBC(3,3) ,CM(3,3) ,GNDB(70) 

DIMENSION X(3) , Y(3) 

C TO CALCULATE THE AREA AND BASIC MATRIX OF THE ELEMENT 

XY - X (2) * Y (3) - Y (2) * X(3) 

YX - X (1) * (Y (2) - Y (3) ) + Y Cl) * (X(3) - X(2)) 

XPY - XY + YX 
EAREA(M) - .5 * ABS (XPY) 

B (M, 1) - Y (2) - Y (3) 

B (M, 2) - Y (3) - Y (1) 

B (M, 3) - Y Cl) - Y (2) 

C (M, 1) - X (3) - X (2) 

C (M, 2) - X (1) - X (3) 

C (M, 3) * X (2) - X Cl) 

ESL(M.I) - SQRT (B (M, 3) **2 + C (M, 3) **2) 

ESL (M, 2) = SQRT (B (M, 1) **2+C (M, 1) **2) 

ESL (M, 3) - SQRT(B(M,2)**2+C(M,2)**2) 

RETURN 
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END 

SUBROUTINE CONCP (X , Y , TEMP , M, DUMK , NLEM) 

Q * it * it it it it it it it it it it it it it it it it it it it * it it it it it it it it it it it it it it it * * it it it it it it it it it it 

C * THIS IS MAIN SUBTROUTINE TO EVALUATE CONDUCTION * 

C * AND CAPACITANCE MATRIX FOR EACH ELEMENT * 

Q it it it it it it it it it it it it it it it it it it it it it it it it it it it it it it it it it it it it it it *********** * * 

COMMON/MATG/ AMC(140,151) ,RM(140) ,THETA(140) ,FLD(140) , 

& AM(140, 140) ,BM(140) ,CIJ(3,3) ,KIJ(3,3) 

& /TEPS/ TEMPR, TEMPI, TMEL,DETP,DELTD, TEMPG1 (140) ,THETAG(140) , 

& THETAFC140) ,THETAO(140) ,THETTO(140) , THETAB (70) , TEMPD 

& /PROC/ HRSL , ALPAR , DIST , B I , KREF , ROUR , SPECR , EMI S , HTC ( 1 40 , 3) , 

£, BOLT, VOLSPR, WIDTH, QE (140, 3) ,NRESIS,RESIS, IX 

& /ELMT/ NEL,NP,GND(140,3) ,B (140,3) , C (140, 3) ,EAREA(140) ,NPP 1 , 

& IES (140 , 2) , ESL (140 , 3) , BETA, XC (140 , 3) , YC ( 140 , 3) , DEFO , 

& NE(140) ,TLENG,PN0D(5,30) ,EBC(3,3) ,CM(3,3) ,GNDB(70) 

DIMENSION X(3) , Y(3) ,TEMP(3) ,KN(3) ,DUMK(140) 

REAL KN 

C TO CALCULATE ELEMENT CONDUCTION M AND CAPACITANCE [C] 

C MATRIX ON ELEMENT BY ELEMENT BASIS 


C TO PROVIDE COMMON MATRIX DATA 

IF (DUMK (M) .EQ. 4.) THEN 
CALL FI CLAY (TEMP ,M) 

ELSE 

C CALCULATE NODAL CONDUCTIVITY AND VOLUMETRIC SPECIFIC HEAT.... 

CALL PROPTY (TEMP, CEE, KN, NLEM, M) 

C CALCULATE CONDUCTION [KIJ] AND CAPACITANCE [CIJ] MATRIES 


CALL CKAN (CEE, KN, TEMP, M) 
END IF 
RETURN 
END 


SUBROUTINE FICLAY (TEMP , M) 


C 

C 

C 

c 

c 

c 


* it it it it it it it it it it it it it it it it it it it it it it it it it it it it it it it * * * it it it it it it it it it it it it it it it it it it 


* THIS SUBROUTINE EVALUATES CONDUCTION AND CAPACITANCE * 

* MATRIX IN FICTITIOUS-LAYER WHICH HAS BEEN ARTIFICIALLY * 

* CREATED TO SMOOTH OUT THE EFFECT OF THE ABRUPT CHANGE * 

* OF TEMPERATURE AT THE BOUNDARY OR FOR HEAT PIPE SHELL. * 


**■ 


’ it it it it it it it it it it it it if 


• * : : it it it it it it it it it it it it it it it it it it it it it it it it it it it it it it it it it it it it it 


IMPLICIT REAL (K) 

COMMON/MATG/ AMC(140, 151) , RM ( 140) , THETA(1 40) ,FLD(140) , 

& AM (140, 140) , BM ( 140) ,CIJ(3,3) ,KIJ(3,3) 

£, /PROC/ HRSL, ALPAR , DIST , BI , KREF , ROUR , SPECR , EMI S , HTC (140, 3) , 

& BOLT, VOLSPR, WIDTH, OE(140, 3) , NRES IS , RES IS , IX 
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& /ELMT/ NEL,NP,GND(140,3) ,B(140,3) ,C(140,3) ,EAREA(140) ,NPP1, 

& IES(140,2) ,ESL(140,3) , BETA,XC (140, 3) , YC (140 , 3) , DEFO, 

L NE(140) ,TLENG,PNOD(5,30) ,EBC(3,3) ,CM(3,3) ,GNDB(70) 

DIMENSION TEMP (3) 

CALL PROSHEL (TEMP , KFIC , CFIC) 

DO 100 I - 1,3 
DO 110 J - 1,3 

EBC(I.J) * (B (M, I) *B (M, J) + C (M, I) *C (M, J) ) /EAREA(M) 

IF (I .EQ. J) THEN 

CM(l, J) - EAREA (M) / 6 . 

ELSE 

CM(I,J) - earea(m) / 12 . 

END IF 

KI J (I , J) - KFIC * EBC(I, J)/4./KREF 
CIJ(I,J) * CFIC * CM (I , J) /VOLSPR 
110 CONTINUE 
100 CONTINUE 
RETURN 
END 

SUBROUTINE PROSHEL (TEMP, KFIC, CFIC) 

q ********************************************** ************ ** 

C * SUBROUTINE CALCULATES CONDUCTIVITY AND SPECIFIC HEAT FOR * 

C * HEAT PIPE SELL BASED ON AVERAGE TEMPERATURE. * 

Q * ** * ** ** * **** ****** * * *************************** ********** ** 


IMPLICIT REAL (K) 

DIMENSION TEMP (3) 

C CONDUCTIVITY AND SPECIFIC HEAT FOR HASTELLAY X. . 

DATA TKl,TK2,TK3,TK4/373. ,573. ,773., 973./ 

& Kl ,K2 ,K3 ,K4 /ll. 1, 14.7,20.6,22.8/ 

& CPl,CP2,CP3/498. ,582. ,699./ 

& TCI, TC2.TC3/588. ,923. ,1143/ 

ROU - 8220. 

TAVG - (TEMP (1) +TEMP (2) +TEMP (3) ) /3 . 

IF (TAVG .LT. TK1) THEN 
KFIC - Kl 

ELSE IF (TAVG .GE. TK1 .AND. TAVG .LT. TK2) THEN 
KFIC - K1+ (K2-K1) * (TAVG-TK1) / (TK2-TK1) 

ELSE IF (TAVG .GE. TK2 .AND. TAVG .LT. TK3) THEN 
KFIC - K2+ (K3-K2) * (TAVG-TK2) / (TK3-TK2) 

ELSE IF (TAVG .GE. TK3 .AND. TAVG .LT. TK4) THEN 
KFIC - K3+ (K4-K3) * (TAVG-TK3) / (TK4-TK3) 

ELSE 

KFIC - K4 
END IF 

IF (TAVG .LT. TCI) THEN 
CFIC - ROU * CPI 

ELSE IF (TAVG .GE. TCI .AND. TAVG .LT. TC2) THEN 
CFIC = ROU * (CP1+ (CP2-CP1) * (TAVG-TC1) / (TC2-TC1) ) 
ELSE IF (TAVG .GE. TC2 .AND. TAVG .LT. TC3) THEN 
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CFIC - ROU * (CP2+ (CP3-CP2) * (TAVG-TC2) / (TC3-TC2) ) 
ELSE 

CFIC » ROU * CP3 
END IF 
RETURN 
END 

SUBROUTINE PROPTY (TEMP , CEE , KN , NLEM , M) 


q is is is it it it is it it it it it it it it it it * it is is is it it it it it it it it it it it it it it it it it it it it it it is it it 

C * THIS SUBROUTINE CALCULATES CONDUCTIVITY AND * 
C * VOLUMETRIC SPECIFIC HEAT AT NODE IN ELEMENT * 
C * WHICH HAS PHASE CHANGE * 

Q it it it it it it it it it it it it it * * * * * it it it it it it * * * * * * * * * * * it it it it * * * * * * it * 


IMPLICIT REAL (K) 

COMMON/TEPS/ TEMPR , TEMP I , TMEL , DETP , DELTD , TEMPG1 (140) , THETAG(IAO) , 
S, THETAF (140) .THETAO (140) , THETTO (140) , THETAB (70) , TEMPO 

S, /FROC/ HRSL,ALPAR,DIST,BI,KREF,R0UR,SPECR,EMIS,HTC(140,3) , 

S, BOLT, VOLSPR, WIDTH, QE(140, 3) ,NRESIS,RESIS, IX 

& /ELMT/ NEL.NP ,GND(140, 3) ,B (140, 3) ,C (140, 3) , EAREA(140) ,NPPl , 

& IES(140,2) ,ESL(140,3) ,BETA,XC(140,3) ,YC(140,3) , DEFO, 

& NE (140) , TLENG , PNOD (5 , 30) ,EBC(3,3) ,CM(3,3) ,GNDB(70) 

& /SPLIN/TIME (3) .TIMER (20) ,QSTAG(20) ,XP1(25) , CPI (4,25) ,ND1, 

L XP2 (60) , CP2 (4 , 60) , ND2 , XP3 (25) , CP 3 (4 , 25) , ND3 , TS , 

£, XP4 (25) ,CP4(4,25) ,ND4,XP5(15) ,CP5(4,15) ,ND5,XP6(15) , 

£, CP6 (4,15) ,XP7 (15) ,CP7 (4, 15) ,ND7,XP8(15) ,CP8(4,15) , 

& ND6 , ND8 ,XP9 (15) ,CP9(4,15) ,ND9,XP0(15) ,CP0(4,15) ,ND0, 

{, ELSS(30) ,TESL(30) 

DIMENSION TEMP (3) ,KN (3) ,KFN (3) ,KLN (3) ,KSN (3) ,ENTP (3) 

TMEP - TMEL + DETP 
TMES - TMEL - DETP 
DO 50 I - 1,3 
KFN(I) - 0. 

KLN(I) - 0. 

KSN(I) - 0. 

50 CONTINUE 

DO 100 I - 1,3 

IF (TEMP (I) .LT. TMES) THEN 
TEMPS-TEMP(I) 

KSN(I) - PCUBIC (TEMPS ,ND5 ,XP5 , CP5) /KREF 
CSN - PCUBIC (TEMPS, ND7.XP7.CP7) 

ROUSN - PCUBIC (TEMPS ,ND6,XP6,CP6) 

ENTP(I) - CSN * ROUSN * TEMP (I) 

ELSE IF (TEMP (I) .GT. TMEP) THEN 
TEMP L= TEMP (I) 

KLN(I) = PCUBIC (TEMPL.ND8 ,XP8 , CP8) /KREF 
CLN - PCUBIC (TEMPL , NDO , XPO , CPO) 

ROULN = PCUBIC (TEMPL, ND9 ,XP9 , CP9) 

CSN = PCUBIC (TMES, ND7,XP7,CP7) 

ROUSN - PCUBIC (TMES ,ND6 , XP6 , CP6) 

ENTP(I) = CSN , ’'ROUSN*TMES J -HRSL 4 -CLN , 'ROULN* (TEMP (I) -TMEP) 

ELSE 

KSS = PCUBIC (TMES, ND5,XP5,C?5) 

CSN = PCUBIC (TMES, ND7,XP7,CP7) 


ORIGINAL PAGE IS 
OF POOR QUALITY 



ORIGINAL PAGE IS 
OF POOR QUALITY 


ROUSN - PCUBIC (TMES ,ND6 ,XP6 , CP6) 

KLP - PCUBIC (TMEP ,ND8 ,XP8 , CP8) 

KFN(I) - (KSS + (KLP - KSS) * (TEMP (I) - TMES) / (2. *DETP) ) /KREF 
ENTP(I) -CSN*ROUSN*TMES + HRSL* (TEMP (I) - TMES) / (2. *DETP) 

END IF 
100 CONTINUE 
DETX - 0. 

DETY - 0. 

DO 200 I - 1,3 

DETX - DETX + (B(M,l)*TEMP(I))/(2.*EAREA(M)) 

DETY - DETY + (C(M, I) *TEMP (I) ) / (2. *EAREA(M) ) 

200 CONTINUE 

IF (DETX .NE. 0. .AND. DETY .NE. 0.) THEN 
ENTPX - 0. 

ENTPY - 0. 

DO 300 I - 1,3 

ENTPX - ENTPX + (B(M,I)*ENTP(l))/(2.*EAREA(M)) 

ENTPY - ENTPY + (C(M,I)*ENTP(I))/(2.*EAREA(M)) 

300 CONTINUE 

IF (NLEM ,EQ. 1) THEN 

CEE - ( (ENTPX**2 + ENTPY**2) / (DETX ,,r *2 + DETY**2) ) /VOLSPR 
ELSE 

CEE - ( (ENTPX*DETX + ENTPY*DETY) / (DETX**2 + DETY**2) ) /VOLSPR 
END IF 
ELSE 

IF (TEMP (1) .LT. TMES) THEN 
CEE - CSN*R0USN/ VOLSPR 
ELSE IF (TEMP (1) .GT. TMEP) THEN 
CEE - CLN ,! ROULN/VOLSPR 
END IF 
END IF 

DO 400 1-1,3 

KN (I) - KFN(I) + KLN(I) + KSN(I) 

400 CONTINUE 
RETURN 
END 


SUBROUTINE CKAN (CEE , KN , TEMP , M) 


C ***************************************** 

C * THIS SUBROUTINE CALCULATES CONDUCTION * 
C * AND CAPACITANCE MATRICES FOR ELEMENT * 
C * WHICH HAS PHASE CHANGE. * 

(2 * * * * * * * * it * * * it it it it it it it ************* * * ****** * 


IMPLICIT REAL (K) 

COMMON/MATG/ AMC(140,151) ,RM(140) , THETA (140) ,FLD(140) , 

& AM (140, 140) , BM (140) ,CIJ(3,3) ,KIJ(3,3) 

& /PROC/ HRSL,ALPAR,DIST,BI,KREF,ROUR,SPECR,EMIS,HTC(140,3) , 

4 BOLT, VOLSPR, WIDTH, QE(140, 3) , NRESIS , RESIS , IX 

& /ELMT/ NEL, NP , GND(140 ,3) ,B(140,3) ,C(140,3) , EAREA (140) ,NPP1 , 

4 IES (140,2) ,ESL(140,3) .BETA, XC ( 140,3) ,YC (140, 3) , DEFO , 

& NE (140) , TLENG, PNOD (5 , 30) ,EBC(3,3) ,CM(3,3) ,GNDB(70) 

4 / SPLIN/TIME (3) , TIMER (20) ,QSTAC(20) ,XP1 (25) ,CP1 (4,25) ,ND1, 

& X?2 (60) , CP2 (4 , 60) ,ND2,XP3(25) ,CP3 (4,25) ,ND3,TS, 
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& XPA (25) , CPA (A, 25) ,NDA,XP5(15) ,CP5(A,15) ,ND5,XP6(15) , 

& CP6 (A , 15) ,XP7 (15) , CP7 (A , 15) , ND7 ,XP,8 (15) ,CP8(A,15) , 

& ND6,ND8,XP9(15) ,CP9(A,15) ,ND9,XP0(15) ,CP0(A,15) ,ND0, 

& ELSS(30) ,TESL(30) 

DIMENSION KN (3) , TEMP (3) 

TAVG - (TEMP (1) + TEMP (2) + TEMP(3))/3. 

VOIDF - .68A 

KSW - PCUBIC (TAVG.ND3 ,XP3 , CP3) 

KSW - KSW/KREF 

CPSW - PCUBIC(TAVG,NDA,XPA,CPA) 

DENSW - 8027. 

CPSW - CPSW*DENSW/VOLSPR 
KEE - 0. 

DO 100 I - 1,3 
KEE - KEE + KN (I) 

100 CONTINUE 

KEE - KEE/3. 

KEE - KEE* ( (KEE+KSW) - (1-VOIDF) * (KEE-KSW) ) / 

& ( (KEE+KSW) + (1-VOIDF)* (KEE-KSW)) 

CEE - (VOIDF*CEE)+ (1-VOIDF) *CPSW 
DO 200 I - 1,3 
DO 210 J - 1,3 

EBC (I , J) - (B(M,I)*B(M,J) + C (M, I) *C (M, J) ) /EAREA(M) 

IF (I .EQ. J) THEN 

CM (I , J) - EAREA (M) / 6 . 

ELSE 

CM (I , J) - EAREA (M) / 1 2 . 

END IF 

CIJ(I,J) » CEE*CM(I,J) 

KIJ(I,J) - KEE*EBC(I, J)/A. 

210 CONTINUE 
200 CONTINUE 
RETURN 
END 


SUBROUTINE THIRDB (THETA C , L , CONF , M) 


C 

C 

C 

C 


it * it it it it it it it it it it it it it it it it it it it it it it it it it it it it it it it it it it it it it it it it it it it it it it it it it it it it it it it it it it it it 


* THIS SUBROUTINE CALCULATES ELEMENT MATRICES DUE TO SURFACE * 

* INTEGRALS FOR CONVECTIVE BOUNDARY CONDITIONS * 


it it it it it it it it it it it it it it it it it it it it it it it it it it it it it it it it it it it it it ' 


■ it it it it it it it it 


' it it it it 


IMPLICIT REAL (K) 

COMMON/MATG/ AMC(1A0,151) ,RM(1A0) ,THETA(1A0) ,FLD(1A0) , 

& AM(IAO.IAO) , BM (1A0) , Cl J (3,3) ,KIJ(3,3) 

& /TEPS/ TEMPR, TEMPI, TMEL t DETP,DELTD,TEMPGl(l AO) .THETAG(IAO) , 

& THETAF (1A0) , THETAO (1A0) , THETTO (1A0) , THETAB (70) , TEMPD 

& /PROC/ HRSL, ALPAR , DIST , BI , KREF , ROUR , SPECR , EMIS ,HTC (1 AO , 3) , 

L BOLT, VOLSPR, WIDTH, QE(1A0, 3) ,NRESIS,RESIS, IX 

& /ELMT/ NEL.NP , GND (1A0, 3) ,B(1A0,3) ,C(1A0,3) ,EAREA(1A0) ,NPP1, 

& IES (1A0, 2) , ESL (1A0 , 3) ,BETA,XC(1A0,3) , YC(1A0,3) , DEFO , 

S, NE(IAO) , TLENG, PNOD (5,30) , EBC (3,3) ,CM(3,3) ,GNDB(70) 

DIMENSION CONF (2 , 3) ,CONK(3,3) , IE(2) 

BI = (DIST/KREF) *HTC (M, IES (M, L) ) 

IF (IES (M, L) .EO. 1) THEN 
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IE(L) - 3 
ELSE 

IE (L) - IES(M,L) - 1 . 

END IF 

DO 100 I - 1,3 

IF (I .EQ. IE (L) ) THEN 
C0NF(L,I) - 0. 

ELSE 

CONF(L.I) - BI * ESL(M,IES(M,L)) * THETAC/2. 
END IF 

DO 110 J - 1,3 

IF (I .EQ. IE (L) .OR. J .EQ. IE(L)) THEN 
CONK (I, J) - 0. 

ELSE IF (I .EQ. J) THEN 

CONK (I , J) * BI*ESL(M,IES(M,L))/3. 

ELSE 

CONK (I, J) - BI*ESL(M,IES(M,L))/6. 

END IF 
110 CONTINUE 

100 CONTINUE 

DO 200 1-1,3 
DO 210 J - 1,3 

KI J (I , J) - KIJ(I.J) + CONK (I, J) 

210 CONTINUE 

200 CONTINUE 
RETURN 
END 


SUBROUTINE SECNDB (L , QEF , M, NFLUX) 


************************************************************** 

* THIS SUBROUTINE CALCULATES ELEMENT MATRICES DUE TO SURFACE * 

* INTEGRALS FOR B.C.'S OF SECOND KIND, AND HEAT FLUX (QE) IS * 

* CONSTANT BETWEEN NODES * 

************************ ************** ************************ 


IMPLICIT REAL (K) 

COMMON/TEPS/ TEMPR , TEMPI , TMEL, DETP ,DELTD, TEMPG1 (140) , THETAG(IAO) , 
& THETAF (140) , THETAO (140) ,THETT0(140) , THETAB (70) , TEMPD 

L /PROC/ HRSL,ALPAR,DIST,BI,KREF,R0UR,SPECR,EMIS,HTC(140,3) , 

& BOLT, VOLSPR, WIDTH, QE(140, 3) ,NRESIS,RESIS, IX 

L /ELMT/ NEL.NP , GND (140, 3) , B (140 , 3) , C (140, 3) , EAREA(140) ,NPP1 , 

5. IES (140,2) , ESL (140 , 3) ,BETA,XC(140,3) , YC(140,3) ,DEFO, 

6, NE (140) ; TLENG , PNOD (5 , 30) ,EBC(3,3) ,CM(3,3) ,GNDB(70) 
DIMENSION IE (2) ,QEF(2,3) 

INTEGER GND 

IF (IES (M, L) .EQ. 1) THEN 
IE (L) - 3 
ELSE 

IE (L) - IES (M, L) - 1 
END IF 


C CONSTANT HEAT FLUX 

C ESL = LENGTH OF SIDE 
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C DLQE - DIMENSIONLESS CONSTANT HEAT FLUX ON SIDE OF ELEMENT 

IF(NFLUX .EQ. 1) THEN 

IF (IES (M, L) .EQ. 1) THEN 
QES - QE (M, IES (M, L) ) 

ELSE 

QS - QE (M, IES (M, L) ) 

CALL FLUXD (M,QES) 

END IF 
ELSE 

QES » QE(M, IES (M, L) ) 

END IF 

DLQE - QES*DIST/ (KREF* (TMEL - TEMPI)) 

IF (IES (M,L) .EQ. 1) THEN 

QEF (L, 1) ” (ESL (M, IES (M, L) ) /2 . ) *DLQE 
QEF(L,2) - QEF(L, 1) 

QEF (L, 3) - 0. 

ELSE IF (IES (M, L) .EQ. 2) THEN 
QEF (L, 1) - 0. 

QEF (L , 2) - (ESL (M, IES (M , L) ) / 2 . )' *DLQE 
QEF(L,3) - QEF(L,2) 

ELSE 

QEF (L, 1) - (ESL (M, IES (M, L) ) /2 . ) *DLQE 
QEF(L,2)«0. 

QEF(L,3) - QEF(L, 1) 

END IF 
RETURN 
END 


SUBROUTINE RADB (L , TS , NTS , RADF , M, RADK) 

Q **************Tflf******************AVc******* , ***********Vc^t****VinV3V* 

C * THIS SUBROUTINE CALCULATES ELEMENT MATRICES DUE TO SURFACE * 

C * INTEGRALS FOR B.C.’S OF RADIATION. * 

Q ************************************************************** 


IMPLICIT REAL (K) 

COMMON/MATG/ AMC(140,151) ,RM(140) ,THETA(140) ,FLD(140) , 

& AM(140, 140) ,BM(140) ,CIJ(3,3) ,KIJ(3,3) 

4 /TEPS/ TEMPR, TEMPI, TMEL, DETP , DELTD, TEMPG1 (140) ,THETAG(140) , 

Sr THETAF (140), THETAO (140), THETTO (140), THETAB (70) , TEMPD 

S, /PROC/ HRSL,ALPAR,DIST,BI,KREF,R0UR,SPECR,EMIS,HTC(140,3) , 

S, BOLT, VOLSPR, WIDTH, QE(140, 3) , NRESIS ,RESIS , IX 

S. /ELMT/ NEL, NP , GND (140 ,3) , B (140 , 3) ,C (140 , 3) , EAREA(140) ,NPP1 , 

S< IES(140,2) ,ESL(140,3) , BETA.XC (140, 3) ,YC(140,3) ,DEFO, 

S, NE (140) , TLENG, PN0D(5 , 30) ,EBC(3,3) ,CM(3,3) ,GNDB(70) 

DIMENSION IE (2) , RADF (2, 3) , RADK (3, 3) 

INTEGER GNDN , GND1 , GND2 , GND , TS 

RADC - (DIST/KREF) *EMIS ,v BOLT 
DELTT = TMEL - TEMPI 
IF (IES (M, L) .EQ. 1) THEN 
IE (L) = 3 
ELSE 

IE (L) = IES (M, L) - 1 
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END IF 

THETAR - (TEMPR - TMEL) / (TMEL - TEMPI) 
C CONSTANT BETA 


ETC - ES L (M , IES (M , L) ) ’'THETAR 
IF (IES (M,L) .EQ. 1) THEN 
GND1 - GND(M, 1) 

GND2 - GND (M,2) 

ELSE IF (IES (M,L) .EQ. 2) THEN 
GND1 - GND (M , 2) 

GND2 - GND(M, 3) 

ELSE 

GND1 - GND(M, 1) 

GND2 - GND (M , 3) 

END IF 

IF (TS .LE. NTS) THEN 

DTHETA - DELTT * (THETAG(GNDl) + THETAG (GND2) ) /2 . 

ELSE 

DTHETA - DELTT * (THETAO (GNDf) + THETAO (GND2) ) /2 . 

END IF 

RBETAC - RADC * ( (DTHETA+TMEL) **2 + TEMPR**2) * (DTHETA+TMEL+TEMPR) 
DO 100 1-1,3 

IF (I .EQ. IE (L) ) THEN 
RADF (L, I) - 0. 

ELSE 

RADF (L, I) - ETC * RBETAC /2. 

END IF 

DO 110 J - 1,3 

IF (I .EQ. IE (L) .OR. J .EQ. IE (L) ) THEN 
RADK(I,J) - 0. 

ELSE IF (I .EQ. J) THEN 

RADK(I,J) - RBETAC*ESL(K,IES(M,L))/3. 

ELSE 

RADK ( I , J ) - RBETAC*ESL(M,IES(M,L))/6. 

END IF 
110 CONTINUE 

100 CONTINUE 

DO 200 I - 1,3 

DO 250 J - 1,3 

KIJ(I,J) - KI J (I , J) + RADK (I , J) 

250 CONTINUE 

200 CONTINUE 
RETURN 
END 


SUBROUTINE FLUXD(M.QES) 


C 

C 

c 

c 


f * * it it ifit't'.tztzt it * it s: r 


it it it it it it it it it it it it it it it it it it it it it it it it it it it it it it it it it it it it it it it it it it it it it it 


* WHEN HEAT FLUX DEPEND ON TIME, CALCULATE HEAT FLUX ON BOUNDARY * 

* SURFACE BY LINEAR INTERRELATION. 


COMMON/ELMT/ NEL, NP , GND (1 40 , 3) ,B (140,3) ,C(140,3) ,EAREA(140) ,NPP1, 
& IES (140, 2) , ESL (140,3) , BETA.XC (140, 3) , YC ( 140, 3) , DEFO , 
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S, NE(IAO) , TLENG , PNOD (5 ,30) ,EBC(3,3) ,CM(3,3) ,GNDB(70) 

L /SPLIN/TIME (3) ,TIMER(20) ,QSTAG(20) ,XP1(25) ,CP1(A,25) ,ND1, 

& XP2 (60) ,CP2(A,60) ,ND2, XP3 (25) , CP3 (A , 25) ,NQ3,TS, 

L XPA (25) ,CPA(A,25) ,NDA,XP5(15) ,CP5(A,15) ,ND5,XP6(15) , 

& CP6 (A, 15) , XP7 (15) , CP7 (A, 15) ,ND7,XP8(15) ,CP8(A,15) , 

& ND6,ND8,XP9(15) ,CP9(A,15) ,ND9,XP0(15) ,CP0(A,15) ,ND0, 

& ELSS (30) ,TESL(30) 

INTEGER NL , NU , NTIME , TS , GND , PNOD , K , KMl , KM2 

IF (TS .LT. 3) THEN 
LTS - TS 
ELSE 
LTS - 3 
END IF 
K - 2 
KMl - 1 
KM2 - KMl - 1 

50 IF (GND(M, 3) .EQ.. .PNOD (1, KMl) .AND. GND (M, 2) .EQ. PNOD(l,K)) 

& THEN 

IF (KMl .EQ. 1) THEN 

XI - TESL (1) / (2 . *TLENG) 

ELSE 

XI - (TESL(KM2) + ELSS (KMl) /2.) /TLENG 
END IF 
ELSE 

K - K + 1 
KMl - K - 1 
KM2 - KMl - 1 
GO TO 50 
END IF 

QSCALE - PCUBIC(XI,ND1 ,XPl ,CP1) 

TMI - TIME (LTS) 

QSFLUX - PCUB IC (TMI , ND2 , XP2 , CP2) 

QES - ABS (QSCALE) *QSFLUX 

RETURN 

END 


SUBROUTINE SPLINE (N, XI, C) 


***************************************************************1' i 

* SUBROUTINE SPLINE USES GAUSS ELIMINATION TO CALCULATE C2,I - * 

* SI WITH GIVEN THE NUMBER S Cl.I-FI, AND C2,I = SI ,C2,N + 1*SN+1* 


DIMENSION XI (60) ,C(A,60) ,D(60) ,DIAG(60) 

DATA DIAG(l) ,D(1)/1.,0./ 

NP - N+l 
DO 10 M - 2 , NP 

D(M) - XI (M) - XI(M-l) 

DIAG(M) - (C(1,M) - C(1,M-1))/D(M) 

10 CONTINUE 

DO 20 M = 2,N 

C(2,M) = 3 . " (D (M) *DIAG (M 4 - 1) + D(M+1) *DIAG (M; ) 
DIAG(M) = 2 . " (D (M) r D(M+1)) 

20 CONTINUE 
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DO 30 M - 2, N 

G - - D(M+1)/DIAG(M-1) 

DIAG(M) - DIAG(M) + G ,C D(M-1) 

C(2,M) - CU.M) + G*C(2,M-1) 

30 CONTINUE 

DO 40 M - 2.N 
NP - NP - 1 

C(2,NP) - (C (2 , NP) -D (NP) *C (2 , NP+ 1) ) /DIAG (NP) 
40 CONTINUE 
RETURN 
END 


SUBROUTINE CALCF(N.XI.C) 


************************************************************ 

* WITH FI STORED IN Cl, I AND SI STORED IN C2,I, SUBROUTINE * 

* CALCULATES C3,I, C4.I ,1 - 1 N * 

************************************************************ 

DIMENSION XI (60), C(4, 60) 

DO 10 I - l.N 

DX * XI(I+1) - XI (I) 

DIVDF1 - (C (1,1+1) - C(1,I))/DX 
DIVDF3 - C(2, 1) + C (2,1+1) - 2.*DIVDF1 
C (3 , 1) - (DIVDF1 - C (2, 1) - DIVDF3) /DX 
C (4, I) - DIVDF3/DX/DX 
10 CONTINUE 
RETURN 
END 


FUNCTION P CUBIC (XBAR , N , XI , C) 

DIMENSION XI (60) ,C(4,60) 

I - 1 

DX - XBAR - XI (I) 

IF (DX) 10,30,20 
10 IF (I .EQ. 1) GO TO 30 
I - I -1 

DX - XBAR - XI (I) ■ 

IF (DX) 10, 30, 30 

19 I - I + 1 
DX - DDX 

20 IF ( I .EQ. N) GO TO 30 
DDX - XBAR - XI(I+1) 

IF (DDX) 30,19,19 

30 P CUBIC - C(1,I)+DX*(C(2,I)+DX*(C(3,I)+DX*C(4,I))) 

RETURN 
END 


SUBROUTINE FORMF (L , CONF , QEF , RADF , M) 


* THIS SUBROUTINE FORMS THE GLOBAL COLUMN VECTORS [FLD] * 


C 

C 

c 
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100 


common/matg/ 

L. 

& /ELMT/ 


& 

& 


AMC (140, 151) ,RM(140) ,THETA(140) ,FLD(140) , 

AM(140, 140) ,BM(140) ,CIJ(3,3) ,KIJ(3,3) 

NEL , NP , GND (140,3) ,B(140,3) ,C (140, 3) , EAREA(140) ,NPP1, 
IES(140,2) ,ESL(140,3) , BETA,XC (140 , 3) , YC(140,3) , DEFO, 
NE (140) , TLENG , PNOD (5 , 30) ,EBC(3,3) ,CM(3,3) ,GNDB(70) 


DIMENSION CONF(2,3) ,QEF(2,3) ,RADF(2,3) 

INTEGER GND 
DO 100 I - 1,3 
II - GND (M , I) 

FLD(II) - FLD(II) + CONF (L, I) + QEF(L.I) + RADF (L, I) 
CONTINUE 
RETURN 
END 


SUBROUTINE FORMGM (P , TS , NTS , M) 

Q * it * * * * * * * * * ******* * * * * * * * * * * * * * * * * * * * * * * * * * 

c * FORMS SYSTEM MATRICES AM(I,J) AMD BM(I) * 

c ********************************* * * * * ** ** * * 


COMMON/MATG/ AMC(140,151) ,RM(140) ,THETA(140) ,FLD(140) , 

& AM (140, 140) , BM (140) ,CIJ(3,3) ,KIJ(3,3) 

& /TEPS/ TEMPR, TEMPI, TMEL,DETP,DELTD,TEMPG1 (140) ,THETAG(140) , 

& THETAF(140) ,THETAO(140) ,THETTO(140) ,THETAB(70) .TEMPD 

& /PRO C/ HRSL, ALPAR , DIST , BI , KREF , ROUR , SPECR , EMIS , HTC (140,3) , 

& BOLT, VOLSPR, WIDTH, QE( 140, 3) ,NRESIS,RESIS , IX 

L /ELMT/ NEL, NP, GND (140, 3) ,B(140,3) ,C(140,3) ,EAREA(140) ,NPP1, 

& IES(140,2) ,ESL (140, 3) , BETA, XC (140,3) ,YC(140,3) , DEFO, 

L NE (140) , TLENG, PNOD (5, 30) ,EBC(3,3) ,CM(3,3) ,GNDB(70) 

DIMENSION THETAIC140) 

INTEGER GND,P,TS 
REAL KIJ 

IF (TS ,LE. NTS) THEN 
DO 100 I - 1,3 
II - GND (M, I) 

THETAI(II) “ (TEMPG1 (II) - TMEL) / (TMEL - TEMPI) 

100 CONTINUE 

DO 200 I -1,3 
II - GND(M, I) 

DO 210 J - 1,3 
JJ - GND(M,J) 

KIJ(I.J) - KIJ (I , J) *WIDTH/DIST 
CIJ(I,J) - CIJ(I, J)*WIDTH/DIST 

Ah(II.JJ) - AM(II.JJ) + (KIJ(I,J)*BETA + CIJ (I , J) /DEFO) 

BM(II) - BM(II) + (CIJ (I , J) /DEFO- (1 .-BETA) * KI J (I , J) ) *THETAI (JJ) 
210 CONTINUE 
200 CONTINUE 

DO 300 I “ 1 , NP 

BM (I) - BM (I) + FLD(I) *WIDTH/DIST 
300 CONTINUE 
ELSE 

DO 400 I = 1,3 
II = GND(M.I) 

DO 410 J = 1,3 
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JJ - gnd(m, j) 

KI J (I , J) - KIJ(I,J)*WIDTH/DIST 
CIJ(I.J) - CIJ(1, J)*WIDTH/DIST 

AM (II , JJ) -AM(II,JJ) + 3 . *KI J (I , J) /A. + Cl J (I , J) /DEFO 
BM (II) - BM (II) + (CIJ(I, J) /DEFO) *THETAO(JJ) 

& - KIJ(I, J)*THETT0(JJ)/4. 

410 CONTINUE 
400 CONTINUE 

DO 420 I - 1,NP 

BM (I) - BM(I) + FLD(l) *WIDTH/DIST 
420 CONTINUE 
END IF 
RETURN 
END 

SUBROUTINE GLOBMAX (TEMP , TEMPC , M , DUMK , NLEM , NFLUX , NRAD , NTS , I FN) 

(;********************Vf***************************5V**************** * 

C* THIS SUBROUTINE ASSUMBLE ELEMENT MATRIX AND BOUNDARY CONDITION * 
C* INTO GLOBAL SYSTEM MATRIX [AMC] EXCEPT SPECIFIED TEMPERATURE. * 

c****************************************************************** 


COMMON/MATG/ AMC(140, 151) ,RM(140) ,THETA(140) ,FLD(140) , 


& AM(140, 140) ,BM(140) ,CIJ(3,3) ,KIJ(3,3) 

& /TEPS/ TEMPR, TEMP I, TMEL,DETP,DELTD,TEMPG 1(140) ,THETAG(140) , 

L THETAF (140) ,THETAO(140) ,THETTO(140) , THETAB (70) , TEMPD 

& /PROC/ HRSL, ALPAR,DIST, BI ,KREF, ROUR, SPECR, EMIS , HTC (140,3) , 

& BOLT, VOLSPR, WIDTH, QE(140, 3) .NRESIS .RESIS , IX 

& /ELMT/ NEL,NP ,GND(140, 3) ,B (140, 3) ,C (140, 3) , EAREA(140) ,NPP1 , 

& IES (140, 2) ,ESL(140,3) ,BETA,XC(140,3) ,YC (140,3) ,DEFO, 

& NE (140) , TLENG,PNOD(5 ,30) ,EBC(3,3) ,CM(3,3) ,GNDB(70) 

&. /SPLIN/TIME(3) .TIMER (20) ,QSTAG(20) ,XP1(25) , CPI (4,25) ,ND1, 

& XP2 (60) , CP2 (4,60) ,ND2,XP3(25) ,CP3(4,25) ,ND3,TS, 

& XP4 (25) , CP4(4,25) ,ND4,XP5(15) ,CP5(4,15) ,ND5,XP6(15) , 

& CP6 (4, 15) ,XP7 (15) ,CP7 (4, 15) ,ND7,XP8(15) ,CP8(4,15) , 

6. ND6 ,ND8 ,XP9 (15) ,CP9(4, 15) ,ND9,XP0(15) ,CP0(4, 15) ,NDO, 

& ELSS (30) , TESL (30) 

INTEGER M , NLEM , L , P , I FN , TS 

DIMENSION TEMP (3) , DUMK (140) , CONK (3, 3) ,RADK(3,3) ,CONF(2,3) , 

& QEF(2,3) ,RADF(2,3) , TEMPC (140 , 3) 

C TO OBTAIN ELEMENT DATA 


CALL C0NCP(X,Y, TEMP, M, DUMK, NLEM) 

£ * * * * * * * * * * * * * * * * * * * ** ****** * * ****** * * * * * * ******** * * * * * * 
C TO EVALUATE SURFACE INTEGRAL (CONVECTION , HEAT FLUX, 

C RADIATION BOUNDARY CONDITIONS) 

Q * * * * * * * * * * * * * * * * * * * * * ******** * * * * * * ******** * * * * * * * * * * * * 


DO 100 I - 1,3 
DO 150 J = 1,3 
CONK (I, J) = 0. 
RADK (I ; J) = 0. 
CONTINUE 


150 
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100 CONTINUE 

DO 200 I - 1,NP 
FLD(I) - -0.0 
200 CONTINUE 

DO 300 L - 1,2 
DO 350 I - 1,3 
CONF (L, I) - 0. 

QEF (L, I) - 0. 

RADF(L.I) - 0. 

350 CONTINUE 

IF(HTC(M, IES(M.L)) .GT. 0.) THEN 

THETAC - (TEMPC (M, IES (M, L) ) - TMEL) / (TMEL - TEMPI) 
CALL THIRDB (THETAC, L,CONF,M) 

END IF 

IF (QE (M, IES (M, L) ) .NE. 0.) THEN 
CALL SECNDB (L , QEF , M , NFLUX) 

END IF 

IF (HTC (M, IES (M, L) ) .GE. 0. .AND. NRAD .EQ. 1) THEN 
CALL RADB (L , TS , NTS , RADF , M , RADK) 

END IF 

C TO ASSEMBLE THE GLOBAL COLUMN VECTOR [FLD] 

CALL FORMFCL, CONF, QEF, RADF, M) 

300 CONTINUE 

C.*...T0 ASSEMBLE ELEMENT MATRIX INTO GLOBAL SYSTEM MATRIX [AMC] 

CALL FORMGM (P , TS , NTS , M) 

RETURN 

END 


SUBROUTINE FIRSTBC (IFN) 

^ *v?**********************yf********vc******vc******** 

C * THIS SUBROUTINE ENTERS THE SPECIFIED BOUNDARY * 

C * TEMPERATURES IN A SYSTEM OF EQUATIONS. * 

q ... v . * * * , r ,v * * * * V: ****** * ****** * * * * * * * * * ****** * * * * * * * * * 

COMMON/MATG/ AMC(140, 151) ,RM(140) ,THETA(140) , FLD (140) , 

4 AM (140, 140) , BM(140) ,CIJ(3,3) ,KIJ(3,3) 

4 /TEPS/ TEMPR, TEMPI, TMEL, DETP.DELTD.TEMPGl (140) ,THETAG(140) , 

4 THETAF (140) ,THETA0(140) ,THETTO(140) ,THETAB (70) , TEMPD 

4 /ELMT/ NEL,NP,GND(140, 3) , B (140, 3) ,C (140,3) ,EAREA(140) ,NPP1 , 

4 IES (140,2) ,ESL(140,3) , BETA,XC (140 , 3) , YC(140,3) , DEFO, 

4 NE(140) , TLENG , PNOD (5,30) ,EBC(3,3) ,CM(3,3) ,GNDB(70) 

INTEGER GNDB 
DO 100 I - 1 , IFN 
II = GNDB (I) 

IF (II .LT. 0) THEN 
GO TO 100 
ELSE 

bm(ii) = THETAB(I) ORIGINAL page l<5 

DO 110 j - 1,NP OF POOR QUALITY 
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IF (II .EQ. J) THEN 
AM (I I , J) - 1. 

ELSE 

BM(J) - BM ( J) - THETAB (I) *AM(J, II) 
AM (II , J) - 0. 

AM(J,II) - 0. 

END IF 
110 CONTINUE 

END IF 
100 CONTINUE 
RETURN 
END 


SUBROUTINE RESIDUL (RMAX) 

£ »V * * * * * * * 1c 1c 1c Id: It 1c 1c * 1c It 1c Iclc 1c 1c Id: * 1c Icicle 1c 1c Icicle 1c 1c * * * Icicle Icicle Icicle 1c 1c 1c 1c 

C * THIS SUBROUTINE CALCULATE RESIDUALS AT EVERY NODES * 

£ * It * * 1: * * * * * * * 1c 1c 1c 1c 1c 1c 1c 1c ********* * * * * * * * * * * ******** * ** ** * * 


COMMON /MATG/ AMC(140,151) ,RM(140) ,THETA(140) ,FLD(140) , 

& AM(140, 140) ,BM(140) ,CIJ(3,3) ,KIJ(3,3) 

& /TEPS/ TEMPR,TEMPI,TMEL t DETP,DELTD,TEMPGl (140) ,THETAG(140) , 

& THETAF (140) ,THETA0(140) ,THETTO(140) .THETAB (70) , TEMPD 

& /ELMT/ NEL.NP ,GND(140, 3) ,B (140, 3) , C (140, 3) ,EAREA(140) ,NPP1 , 

& IES (140,2) ,ESL(140,3) ,BETA,XC(140,3) , YC(140,3) , DEFO, 

& NE (140) , TLENG , PNOD (5,30) ,EBC(3,3) ,CM(3 , 3) , GNDB (70) 

DO 100 I - 1 ,NP 
DO 110 J - 1 ,NP 

RM(I) - RM(I) + AM (I, J)*THETA(J) 

110 CONTINUE 

RM (I) - RM(I) - BM(I) 

100 CONTINUE 

RMIN - 10000. 

RMAX - 0. 

DO 200 I - 1,NP 

IF(ABS(RM(I)) .GT. RMAX) THEN 
RMAX - ABS(RM(I)) 

END IF 

IF(ABS(RM(I)) ,LT. RMIN) THEN 
RMIN - ABS (RM (I) ) 

END IF 
200 CONTINUE 
RETURN 
END 


SUBROUTINE CHLSKY (TS , NTS , THETAX) 


£ y r y. y- y f y- y- y. y- y f y- y. y. ju y^ y, y. y ; y. y. y. y- y f y. y. y ? y. y ? y ; y f y. y, y. y f y. y- y. y. y. y. y- y. y. y. y- y. yu y. 

C * THIS SUBROUTINE SOLVES THE SYSTEM OF EQUATION * 
C * BY USING CHOLESKY DECOMPOSITION METHOD * 


COMMON /MATG/ AMC(140, 151) ,RM(140) ,THETa(140) ,FLD' 140) , 
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& 

L 

L 

& 

L 

& 


AM(140, 140) ,BM(140) ,CIJ(3,3) ,KIJ(3,3) 

/TEPS/ TEMPR , TEMPI , TMEL , DETP , DELTD , TEMPG1 (140) , THETAG (140) , 
THETAF (140) , THETAO (140) , THETTO (140) , THETAB (70) , TEMPO 
/ELMT/ NEL , NP , GND (140,3) ,B(140,3) ,C(140,3) ,EAREA(140) ,NPPl, 
IES (140,2) ,ESL(140,3) , BETA.XC (140 , 3) , YC (140, 3) , DEFO, 
NE (140) , TLENG , PNOD (5,30) ,EBC(3,3) ,CM(3,3) ,GNDB(70) 


110 

100 


DIMENSION THETAX (140) 
INTEGER TS 
DO 100 I - 1,NP 
DO 110 J - 1 ,NP 
AMC ( I , J ) - AM(I.J) 
CONTINUE 
CONTINUE 


C. .TO SUBSTITUE RM(I) OR BM(l) INTO THE NPP1TH COLUMN OF MATRIX AMC(NP.NPPl) , . 


DO 120 I “ 1,NP 
.IF (TS .LE. NTS) THEN 

AMC(I,NPP1) - AMC(I.NPPl) - RM(I) 

ELSE 

AMC (I ,NPP1) - AMC(I.NPPl) + BM(l) 

END IF 
120 CONTINUE 

C .... TO CALCULATE FIRST ROW OF UPPER UNIT TRIANGULAR MATRIX.... 


DO 200 J - 2,NPP1 

AMC ( 1 , J ) - AMC (1,J) /AMC (1,1) 

200 CONTINUE 

C TO CALCULATE OTHER ELEMENTS OF U AND L MATRICES 


DO 300 I - 2,NP 
J - I 

DO 310 II - J ,NP 



SUM - 0. 

JM1 - J - 1 
DO 320 K - 1 , JM1 

SUM - SUM + AMC (II, K) 

* AMC(K.J) 

320 

CONTINUE 



AMC (II, J) - AMC (II, J) - 

SUM 

310 

CONTINUE 



IP1 - I +‘l 
DO 330 JJ - IP1.NPP1 
SUM =* 0. 

IM1 - I -1 
DO 340 K = 1 , IM1 

SUM - SUM + AMC (I ,K) *AMC (K, J J) 

340 CONTINUE 

AMC (I , J J) = (AMC(I.JJ) ~ SUM) /AMC (I, I) 

330 CONTINUE 
300 CONTINUE 

C TO SOLVE FOR THETAX (I) BY BACK SUBSTITUTION ORIGINAL PAGE IS 

OF POOR QUALITY 
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THETAX(NP) - AMC(NP.NPPl) 

L - NP - 1 
DO 400 NN - 1,L 
SUM - 0. 

I - NP - NN 
IP1 - I + 1 
DO 410 J - IP1, NP 

SUM - SUM + AMC(I,J) * THETAX(J) 
410 CONTINUE 

THETAX(I) - AMC(I.NPPl) - SUM 
400 CONTINUE 
RETURN 
END 


SUBROUTINE INTFLUX ( JFNT , JNT , TVU , TEMPV , TEMPS , QEI) 


Q it * * * * * it it * * it it it it * it it it it it it it it it it it it it it it it it it it it it it it it it it it it it it it it it it it it it * * it it it it it it it it 

C * THIS SUBROUTINE CALCULATE HEAT FLUX AT LIQUID - VAPOR INT- * 
C * FACE BY USING KINETIC THEORY WITH VAPOR AND INTERFACE TEMP* 

Q ************************************************************** 


DIMENSION GNDI (70) , TEMPS (70) .TEMPV (70) , QEI (70) ,MASW(70) 
INTEGER GNDI 
REAL MASW 

KP - 0 

DO 50 I - 1 , JFNT 
QEI (I) - 0. 

. 50 CONTINUE 

DO 100 I - 1 , JNT 
TW - TEMPS (I) 

TV - TVU 

CALL NAVPROP (TW , DENV , VISV , HFG , HG , PSATI , KP) 

CALL NAVPROP (TV , DENV , VISV , HFGV , HG , PSATV , KP) 

MASW (I) - (2.*.7*(23./(2.*3.1416*8314.))**.5)*(PSATI/TW**.5 
& - PSATV/TV**.5) 

QEI (I) - MASW (I) * HFG 
100 CONTINUE 
RETURN 
END 


SUBROUTINE COUPLE ( JNT, QT, TVU, JFNT) 


q it it * * * * * * * * * * ****** * * * * * * * * * ******** * * * * * * ****** * * * * * * * * * * * * * * * 

C* THIS SUBROUTINE CALCULATES THERMAL RESISTANCE IN THE VAPOR * 
C* SPACE BY USING KNOWN HEAT FLUX AND EVALUATES NEW HEAT FLUX * 
C* AND VAPOR TEMPERATURE * 

^ ^ ^ ^ y* y* y* y* y* y* y* y* y« y» y* y* y# y# y^ y* y- y* y* y* y« y» y* y» y# y» y# y# y* y* y« y* y» y^ y# y^ y* y« y« y« y» y« y» y« y* y« y# y# 


COMMON /MATG/ AMC (140, 15 1) ,RM( 140) .THETA (140) ,FLD( 140) , 

5. AM(140, 140) ,BM(140) ,CIJ(3,3) ,KIJ(3,3) 

6. /TEPS/ TEMPR , TEMPI , TMEL, DETP , DELTD, TEMPO 1 (140) , THETAG(140) , 

& THETAF(140) ,THETA0(140) .THETTO(140) , THE TAB (70) .TEMPO 



n n o n o o 
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& 

L 

St 

St 

& 

L 

& 

& 

& 

& 

& 

& 

4 

& 

4 

& 


/PROC/ HRSL, ALPAR , DIST , BI , KREF , ROUR, SPECR , EMIS ,HTC (140,3) , 
BOLT , VOLSPR , WIDTH , QE ( 140 , 3) , NRES I S , RES I S , IX 
/ELMT/ NEL , NP , GND (140,3) ,B(140,3) ,C(140,3) ,EAREA(140) ,NPP1, 
IES (140 , 2) ,ESL(140,3) ,BETA,XC(140,3) ,YC(140,3) ,DEFO, 
NE (140) , TLENG , PNOD (5,30) ,EBC(3,3) ,CM(3,3) ,GNDB(70) 
/INOUT/TEMPS (70) ,TEMPV(70) ,TEMPVS(70) ,QEI(70) ,VX(70) , 

GNDI (70) , VY1 (70) ,VY2(70) ,VY3(70) ,ELSI(30) , 

VY4 (70) , VY5 (70) ,VY6(70) ,HFG,TEMPF (70) 

/SPLIN/TIME (3) .TIMER (20) ,QSTAG(20) ,XP1(25) , CP 1(4,25) ,ND1, 
XP2 (60) , CP2(4,60) ,ND2,XP3(25) ,CP3(4,25) .ND3.TS, 

XP4 (25) ,CP4(4,25) ,ND4,XP5(15) ,CP5(4,15) ,ND5,XP6(15) , 
CP6 (4,15) ,XP7 (15) ,CP7 (4, 15) ,ND7,XP8(15) ,CP8(4,15) , 
ND6,ND8,XP9(15) ,CP9(4,15) ,ND9,XP0(15) ,CP0(4,15) ,ND0, 
ELSS (30) ,TESL(30) 

/OUTPT/MO , IDAY , I YEAR , TITLE , CASE , NUMBER , NO , TEMPCC , HTCC , IOUT , 
TIMEN , DELT , DELT1 , DELTP 
INTEGER JNT , IX , JNTL , JFNT 
REAL JCOB 
DIMENSION QT (70) 


JNT - NUMBER OF NODE WHOSE TEMPERATURE IS GREATER THAN TSTAR 
JNTL - NUMBER OF NODE WHOSE TEMPERATURE IS GREATER THAN TVU 
TEMPD - TEMP. DROP IN VAPOR SPACE [K] 

RESIS “ RESISTANCE IN VAPOR SPACE [K/W] 

RESISN - ARTIFICIALLY ADDED RESISTANCE AT EACH NODE 
TEMPF(I) = TEMP. DROP DUE TO NEW RESISTANCE AT INTERFACE 


QTAL =0.0 

IF (NRESIS .EQ. 2) THEN 
DO 100 I = 1 , JNT 
QEI(I) » QT (I) 

100 CONTINUE 
END IF 
NRESIS = 2 
DO 150 I = 1 , JNT 

IF (TEMPS (I) .GT. TVU) THEN 
JNTL = I 

ELSE IF (TEMPS (I) .LT. TVU .AND. TEMPS(I-I) .GT. TVU) 
& THEN 

XII » (TEMPS (JNTL) - TVU)*ELSI(I) 

XII = XII/ (TEMPS (JNTL) -TEMPS (I)) 

END IF 
150 CONTINUE 

DO 200 I = 1 , JNT 

IF(QEI(I) .GT. 0.) THEN 

QTAL = QTAL * QEI (I) *ELSI (I) *WIDTH 
END IF 

200 CONTINUE 

QEI (JFNT) = 0. 

IF (IX .LE. 1) THEN 

CALL HPVAPOR (JNT, TS , DIST, IND, TIME) 

TEMPD = VY4 (1) - VY4 (JNT) 

RESIS = TEMPD/QTAL 
RESIS = RESIS/23 . 

END IF 


ORIGINAL PAGE IS 
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IF (IX .IE. 100) THEN 

IF (TEMPD .LE. 20.) THEN 
RESISN - RESIS 
ELSE 

RESISN - RES IS* IX*. 01 
END IF 
IX - IX + 1 
IF (IX .EQ. 100) THEN 
IX - 1 
END IF 
END IF 

DO 250 I - 1, JNT 

TEMPF(l) - RESISN*QEI (I) *WIDTH*ELSI (I) 

TEMPS (I) - TEMPS (I) - TEMPF(I) 

250 CONTINUE 

500 RESD - 0. 

TELSI - 0. 

RESDC - 4.5E6*1.A*.021*2.29E11 
DO 260 I - 1 , JNT 

RESDI - ELSI (I) *10. **(-5567. /TEMPS (I)) /TEMPS (I) 
RESDV - ELSI (I) *10.** (-5567. /TVU) /TVU 
RESD - RESD + RESDC* (RESDI - RESDV) 

TELSI - TELSI + ELSI (I) 

260 CONTINUE 

JCOB - (1 - 12818. 5/TVU)*TELSI 
JCOB - JCOB * 10. **(-5567. /TVU) 

JCOB - RESDC* JCOB/ (TVU**2) . 

TVO - TVU 

TVU - TVO - RESD/ JCOB 
IF (ABS (RESD) .LT. .01) THEN 

CALL INTFLUX ( JFNT, JNT , TVU , TEMPV , TEMPS , QEI) 

DO 270 I - 1 , JNT 
QT (I) - QEI (I) 

270 CONTINUE 

ELSE 

GO TO 500 
END IF 
RETURN 
END 


FOLLOWING SUBROUTINES ARE USED TO CALCULATE PRESSURE, VELOCITY, 
DENSITY, TEMPERATURE, AND QUALITY IN VAPOR SPACE 

SUBROUTINE HPVAPOR ( JNTD, TS , DDIST , IND, TIME) 


C ************************************************************************* 

C * THIS PROGRAM CALCULATES THE VAPOR TEMPERATURES, PRESSURES, DENSITIES, * 
C * QUALITY, AND VELOCITY WITH CONSIDERATION OF COMPRESSIBILITY BY USING * 
C * THE RUNGE-KUTTA METHOD. 


C 


variables 



n n n n n 
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C • 

C DIST HEIGHT OF THE VAPOR SPACE [M] 

C INTV INCREMENT IN AXIAL DIRECTION [M] 

C HASH MASS FLOW AT WALL [KG/M**2] 

C TW WALL TEMPERATURE [K] 

C JNT NUMBER OF POINTS WHERE DATA ARE CALCULATED 

C QW HEAT FLUX [W/M**2] 

COMMON/VAPOR/ REYW (70) , II , JNT ,DIST, TW , QW ,MASW, PSAT, YIN, Y2N , 

& Y3N.Y4N, Y6N 

& /INOUT/ TEMPS (70) , TEMPV (70) , TEMPVS (70) , QEI (70) , VX (70) , 

£, GNDI (70) , VY1 (70) , VY2 (70) ,VY3(70) ,ELSI(30) , 

& VY4 (70) , VY5 (70) ,VY6(70) ,HFG, TEMPF (70) 

DIMENSION C (24) ,W(5,9) ,TWI (70) ,Y (6) , MACH (70) ,TIME(3) 

INTEGER II ,N,NW, IND, GNDI , TS 

REAL MASW.MASWN, INTV, XEND, MACH, JCOB 

EXTERNAL FCNl 


C INITIAL PARAMETERS ; 

N - NUMBER OF DIFFERENTIAL EQUATIONS 
Y (1) - PRESSURE OF VAPOR [N/M**2] 

Y (2) - VELOCITY [M/SEC.] 

Y (3) - DENSITY [KG/M**3] 

Y (4) - TEMPERATURE [K] 

C Y(5) - QUALITY OF VAPOR 

C Y (6) = SPECIFIC VOLUME [M**3/KG] 


DIST - DDIST 
JNT - JNTD 
PH “ 0. 

N - 5 
NW - N 
IND - 1 
TOL - .0001 

C CALAULATE INITIAL VALUES 

CALL INITV 
Y (4) - TEMPV (2) 

Y (1) - 2. 29E11/ (Y (4) **.5) *10.** (-5567 . /Y (4) ) 

Y (3) - 23 . *Y (1) / (8314. *Y (4) ) 

Y (5) - 1. 

Y (6) - l./Y (3) 

KP - 1 

CALL NAVPROP (TEMPV (2) , DENV , VISV ,HFG, HG, PSATI , KP) 

Y (2) - QEI(1)*ELSI(1)/(HFG*DIST*Y(3)) 

C SET REFERENCE, VALUE TO NORMALIZE 

YIN = Y ( 1) 

Y2N » Y (2) 

Y3N = Y (3) 

Y4N = Y (4) 
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Y6N » Y (6) 

VX(1) = 0.0 
VX(2) = ELSI(l) 
VY2 (1) - 0.0 
VY2 (2) - Y (2) 

DO 100 I - 1,2 
VY1(I) - Y (1) 
VY3 (I) - Y (3) 
VY4 (I) - Y (4) 
VY5 (I) - Y (5) 
VY6 (i) - Y (6) 
100 CONTINUE 

TWI(l) - TEMPS (1) 
TWI (2) - TEMPS (2) 


C CALCULATE P,V, DENSITY, T, QUALITY BY USING DVERK 

DO 200 II - 3,JNT + 1 
IM1 - II - 1 
QW - QEI(IMI) 

5 PH - PH + 1 


IF (PH .NE. 1 . ) THEN 
IND - 1 
END IF 

INTV - ELSI(IMl) 

XEND - VX(II-l) + INTV 
VX(II) - XEND 
X = VX(IMl) 

C NORMALIZE THE DEPENDANT VARIABLES 

Y (1) - VY1 (IMI) /YIN 

Y (2) - VY2(IM1)/Y2N 

Y (3) - VY3(IM1)/Y3N 
Y (4) - VY4 (IMI) /Y4N 
Y (5) - VY5 (IMI) 

Y (6) - VY6 (IMI) /Y6N 
IF (PH .EQ. 1.) THEN 
TW - TEMPS (II) 

ELSE 

TW - TWI (II) 

END IF 

CALL DVERK (N, FCNL.X, Y , XEND , TOL , IND , C , NW , W , IER) 

IF (IND .LT. 0 .OR. IER .GT. 0) THEN 
WRITE (4,*) IND, IER 
GO TO 300 
END IF 

Y1 - Y(1)*Y1N 
Y2 - Y (2) *Y2N 
Y3 - Y (3) *Y3N 
Y4 = Y (4) *Y4N 

MASWN - 2 . 7* (23 . / (2 . *3 . 1416*83 14 . ) ) ** . 5* (PSAT/TW** . 5 
& - Yl/Y4**.5) 

DMAS = (MASW - MASWN) / (MASW) 

IF (ABS (DMAS) .LT. .02) THEN 
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VYI(II) - Y (1) *Y1N 
VY2 (II) - Y (2) *Y2N 
VY3 (II) - Y (3) *Y3N 
VY4 (II) - Y (4) *Y4N 
VY5 (II) - Y(5) 

VY6 (II) - 1./VY3(II) 

TEMPV(II) - VY4 (II) 

TWI (II) - TW 
PH » 0. 

ELSE 

HFG - 4636437. - 180.82*TW 
RESD - MASW - MASWN ' 

JCOB - 180.82*QW/(HFG**2) 

JCOB - JCOB - 2. *.7* ((23./ (2. *3. 1416*8314.)) **.5) 
i, * (2.29E11* (12818 .8/TW - 1 .) * (10.** (-5567 . /TW) ) /TW**2) 

TW - TW - (RESD/ JCOB) 

TWI (II) - TW 
GO TO 5 
END IF 
200 CONTINUE 
300 CONTINUE 

IJNT - 2 ,,t JNT/3 
WRITE (6, 400) TIME (3) 

WRITE (6 , 410) QW.DIST, INTV ,REYW (2) ,REYW(IJNT) 

WRITE (6, 420) 

WRITE (6, 430) 

DO 350 II - l.JNT+l 

MACH (II) - VY2(II)/SQRT(1. 4*8314. *VY4 (II) /23 . ) 

WRITE (6 , 440) VX(II) , VY1 (II) ,VY2(II) ,VY3(II) ,VY4(II) , 

& VY5(II) ,MACH(II) 

350 CONTINUE 

400 FORMAT (/2X, 83 ('-’) /28X, ’**TIME « ' ,F9.3,2X, ’SEC**'/) 

410 FORMAT (4X, ’HEAT INPUT - ’ , F8 . 1 ,2X, 'DIST - ' ,F6.4,2X, ' INTV - 
& F5. 3, 3X, 'REYNOLD - ' ,F6. 1.2X.F6. 1/) 

420 FORMAT (8X, ’XL (I) ' ,6X, ’PRESSURE' ,4X, 'VELOCITY' ,5X, 

£, 'DENSITY' ,5X, 'TVAP' ,6X, 'QUAL' ,5X, 'MACH'/) 

430 FORMAT (10X, 'M' ,9X, 'N/M**2 ' , 7X, 'M/SEC' ,6X, ' KG/M**3 ’ , 7X, 'K' , 

S. 30X./2X, 83 ('-')/) 

440 FORMAT (5X,E9 . 2, 2X,E11 .5 , 2X,E10. 3 , 3X,E10. 4, 2X, F7. 2, 3X,F6 . 3 , 3X, F6 . 3) 
RETURN 
END 


SUBROUTINE INITV 

q y c * * * * * if * * * * * * * * * * * * * * is is is is it is is is is it it it it it it ******* * * * * * * * * * * * * * 

C * THIS SUBROUTINE CALCULATES INITIAL VAPOR TEMPERARURE. * 

£ A * * * is it is it it it is * * * ****** * is * * * * * * * is * * * * * ********** * * * * * * * * * 

COMMON/ INOUT/ TEMPS (70) , TEMPV (70) ,TEMPVS(70) , QEI (70) , VX (70) , 
GNDI (70) , VY1 (70) , VY2 (70) , VY3 (70) , ELSI (30) , 

& VY4 (70) , VY5 (70) ,VY6(70) , HFG , TEMPF (70) 


INTEGER GNDI 
REAL MASW, JCOB 
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KP - 0 
PHI - 3.1416 
UGAS - 8314. 

TWI - TEMPS (2) 

TV I - TWI - .1 

CALL NAVPROP (TWI , DENV , VIS V , HFG , HG , PSATW , KP) 

100 CALL NAVPROP (TVI,DENV,VISV,HFGV,HG,PSATV,KP) 

NN - NN + 1 

MASW - 2.*.7*(23./(2.*PHI*UGAS))**.5*(PSATW/TWI**.5 
& - PSATV/TVI**.5) 

QEIN - MASW*HFG 
DQEI - (QEIN - QEI (1) ) /QEIN 
IF (NN .LT. 5) THEN 
IF (ABS (DQEI) .LT. .01) THEN 
TEMPV(l) - TV I 
TEMPV (2) - TV I 
ELSE 

RESD - QEI (1) - QEIN 

JCOB - 2. *.7*HFG* (23. / (2. *PHI*UGAS) ) **.5 

JCOB - JCOB * 2.29Elin0.**(-5567./TVI)/(TVI**2) 

JCOB - JCOB*(12818.5/TVI-l.) 

TV I - TV I - RESD/ JCOB 
GO TO 100 
END IF 
END IF 
RETURN 
END 


SUBROUTINE FCN1 (N,X, Y, YPR) 

Q ************************************************************* 

C * THIS SUBROUTINE PROVIDE DERIVITIVE OF GOVERNING EQUATIONS * 
£ ****** ************************ ******************************* 


COMMON/VAPOR/ REYW(70) , II , JNT ,DIST, TW , QW.MASW ,PSAT, YIN , Y2N , 

& Y3N, Y4N, Y6N 

& / INOUT/ TEMPS (70) , TEMPV (70) , TEMPVS (70) , QEI (70) , VX (70) , 

& GNDI (70) , VY1 (70) , VY2 (70) , VY3 (70) ,ELSI(30) , 

& VY4 (70) , VY5 (70) ,VY6(70) ,HFG, TEMPF (70) 

DIMENSION YPR (5), Y (6) 

REAL MASW 

IM1 - II - 1 
Y (6) = VY6 (IMl) 

KP - 0 

CALL NAVPROP (TW , DENW , VISW , HFG , HGW , PSAT , KP) 

MASW - QW/HFG 
VW = MASW/ (DENW) 

REYW(IMl) - - DENW*VW*DIST/VISW 
HCP = 904. 

C THIS SUBROUTINE CALCULATES ALP AR , BETA, F2P0, AND F2P1 


CALL FACTOR (ALPAR , BETA, F2P0, F2P1) 
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THIS SUBROUTINE CALCULATES THE PROPER HIES OF SODIUM, 

TEMPK - Y(4)*Y4N 
KP - 1 

CALL NAVPROP (TEMPK , DDEN , V I S V , HFGV , HG , DP AST , KP) 

Y1 - Y (1) *Y1N 
Y2 - Y(2)*Y2N 
Y3 - Y (3) *Y3N 
Y4 - Y (4) *Y4N 
Y6 - Y (6) 

RSTA - 8314. /23. 

SPEV - 1./Y3 

SPEVF - l./((.927 - . 238E-3* (Y4 - 373))*1000.) 

IF (Y (5) .LE. 0.) THEN 
SPEVG - 0. 

ELSE IF (Y (5) .GT. 1.) THEN 
Y (5) - 1. 

SPEVG - SPEV 
SPEVF - 0. 

ELSE 

SPEVG - (SPEV - SPEVF) /Y (5) + SPEVF 
END IF 

SPERI - (SPEVG - SPEVF) /SPEV 
RTH - (1. - RSTA*Y4/HFGV) 

FIC - 8 . *VISV* (F2P0 - F2P1)/(Y3*DIST*Y2) 

PY1 - ALPAR*MASW/ (BETA*Y2*DIST) 

PY2 - 2. *HFGV' 

PY3 - SPERI* (HGW-HG+(BETA*Y2**2)/2.+(VW**2)/2.) 

PY4 - HFGV/ (BETA*SPEV) 

PY5 - SPERI / SPEV* Y2**2 
PY6 - - PY1* (PY2 + PY3) 

PY7 - PY6 - (PY4+PY5) *FIC/ (8.*DIST) 

PY8 - ALPAR*HFGV*Y(5)*SPEVG*RTH/(BETA*Y1*SPEV**2) 

PY9 - ALPAR*SPERI*HCP*RSTA*Y4**2/(BETA*SPEV*HFGV*Y1) 
PY10 - SPERI + HFGV/ (BETA* Y2**2) - PY8 - PY9 
PY - PY7/PY10 
YPR(l) ** PY/Y1N 

.CALCULATE THE DERIVITIVE OF QUALITY 

QUALl * (Y2**2) *Y (5) *SPEVG*RTH/ (Yl*SPEV**2) 

QUAL2 - (-1./ALPAR + QUALl) *PY 
QUAL3 - FIC*Y2**2/(8.*DIST*SPEV*ALPAR) 

QUAL4 - SPERI*Y2**2/SPEV 

QUAL5 - (QUAL2 - QUAL3)/QUAL4 

QUAL6 - 2.*MASW*SPEV/(Y2*DIST*SPERI) 

QUAL - QUAL5 - QUAL6 
YPR (5) = QUAL 


CALCULATE THE DERIVITIVE OF VELOCITY 

VEL1 = MASW*SPEV/DIST 
VEL2 = Y2*SPERI*QUAL 

VEL3 = Y2*Y(5)*SPEVG*(-RTH)*PY/(Y1*SPEV) 
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VEL - VELl + VEL2 + VEL3 
YPR (2) - VEL/Y2N 

C CALCULATE THE DERIVITIVE OF DENSITY.. 

ROU1 - SPERI*QUAL/ SPEV 

ROU2 - SPEVG*Y(5)*RTH*PY/(Y1*SPEV**2) 

ROU - -ROU1 + ROU2 

YPR (3) - ROU/Y3N 

TY - PY*RSTA*Y4**2/(Y1*HFGV) 

YPR (4) - TY/Y4N 

RETURN 

END 


SUBROUTINE NAVPROP (TEMPK, DENY , VISV , HFG , HG , PSAT , KP) 


Q it * it it it it ******** * * * * * * * * ** *********** * * * * * * * * * * * ** * * * * * * * * * * * * * * * * * * ** 

C * THIS SUBROUTINE CALCULATES THE PROPERITIES OF THE SODIUM SUCH AS * 
C * THE SATURATION PRESSURE, THE DENSITY, THE VISCOSITY, THE LATENT * 
C * HEAT OF VAPORIZATION, AND THE ENTHALPY. * 

Q ******************************************************************** 


T - TEMPK 

IF (KP .EQ. 0) THEN 

PSAT - (2.29Ell/T**.5)*10.**(-5567./T) 
ELSE 

PSAT - 0. 

END IF 

DENV - 2.766E-3*PSAT/T 

VISV - 6.083E-9*T + 1.2606E-5 

HF - 98.973 + 1.4367 * (T - 273.16) 

HF - HF - 2.902E-4* (T - 273.16)**2 
HF - 1000.*(HF + 2. 4E4*EXP (-1 .36E4/T) ) 
HFG - 182. * (25474.93 - .9935*T) 

HG - HF + HFG 

RETURN 

END 


SUBROUTINE FACTOR (ALPAR , BETA, F2P0, F2P1) 


* * * * * * * it it it it it it it it it it it it it it it it it it it * * it it it it it i 


V * V it it it it is ic it i% is i% it i% it 


C * THIS SUBROUTINE PROVIDES THE VALUES OF ALPAR , BETA, F2P0 , * 
C * AND F2P1 CORRESPONDING TO THE GIVEN REYNOLDS NUMBERS. * 

£ it it it it it it it it it * * it it it it it itit it it it it it it it it it it it it it it it it it it it it it it it it it it it it it it it it it it it it it it it it it 


COMMON/ VAPOR/ REYW (70) , I I , JNT , DIST , TW , QW , MASW , PSAT , Y IN , Y2N , 
& Y3N, Y4N, Y6N 

REYW1 = REYW(II-l) 

REYW2 = REYW1 — 2 

IF (REYW1 .LT. 2. .AND. REYW1 . GT. -30.) THEN 

F2P0 - 6.0995 - ,42198*REYW1 - 3.8013E-3 * REYW2 
F2P 1 = - 5.6405 - . 2349 "REYW 1 - 5 . 2913E-3*REYW2 
ALPAR = 1.2049 - 1.0386E-3*REYW1 
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BETA - 1.5574 - 3. 1837E-3*REYW1 
ELSE IF(REYW1 .GE. 2.) THEN 

F2P0 - 6.0995 - .42198*REYW1 - 3.8013E-3 * REYW2 
F2P1 - -6.548 - . 1417*REYWl - . 1233*REYW2 
ALPAR - 1.22 - 1.5082E-2 * REYWl + 2.6689E-3*REYW2 
BETA - 1.612 - 5.0904E-2 * REYWl + 8.6570E-3*REYW2 
ELSE IF (REYWl .LE. -30.) THEN 
F2P0 - 15.34 - .223* (REYWl + 30) 

F2P1 - - 3.06 
ALPAR - 1.227 
BETA - 1.63 
END IF 
RETURN 
END 


SUBROUTINE NEWRA ( JNT , KN , QEI , QEIN, TEMPS , TEMPV , QW) 

Q ******************************* ** ************** ** * * * * * * * * * 

C * CALCULATE NEW WALL TEMPERATURE BY USING NEWTON-RAPHSON * 
C * METHOD FOR NEXT ITERATION.. * 

Q ********************************************************** 


DIMENSION TEMPS (70) , TEMPV (70) , QEI (70) , QEIN (70) 

REAL MASW , MAS WN , JCOB 

DO 100 I - 2, JNT+1 
TW - TEMPS (I) 

TV - TEMPV (I) 

MASW - QEI (I) 

MASWN - QEIN (I) 

RESD - MASW - MASWN 

JCOB - -364. *.7*(SQRT(23./(2. *3. 1416*8314.))) 

& * ( (2 . 29E1 1* (10 . ** (-5567 . /TW) ) /TW) 

& *(-.9935 + (25474.93 - .9935*TW) *( (12818 . 5-TW) /TW**2) ) 

& + .9935* (2. 29E11/TV) *10. ** (-5567. /TV) ) 

TW - TW - RESD/ JCOB 
100 CONTINUE 
RETURN 
END 


SUBROUTINE DATAIN (TA, JK, JNT , TIMEN , 
£, FIRTS , TVU , QT , TEMPD) 


Q * * * * * * ***** *** * ***** * * * * ** * * * * * * * * * * ** *** * * * ** * * * * * * * * * * * * * * * * * * 

C* THIS SUBROUTINE READS IN INTERMEDIATE DATA FOR RESTARTING OF * 
C* MAIN PROGRAM FROM LAST RUN * 

Q ******** ** * * ************************ * ******** *** * * * * * * ****** * * * * 


COMMON/PROC/ HRSL,ALPAR,DIST,BI,KREF,ROUR,SPECR,EMIS,HTC(140,3) , 

L BOLT, VOLSPR, WIDTH, QE(140, 3) , NRESIS , RESIS , IX 

i /ELMT/ NEL,NP ,GND(140, 3) ,B(140,3) , C ( 140 , 3) , EAREA (140) ,NP?1, 

L IES (140,2) , ESL (140 , 3) , BETA, XC ( 140 , 3) ,YC(140,3) , DEFO , 

L NE (140) , TLENG, PNOD (5 ,30) ,EBC(3,3) ,CM(3,3) ,GNDB(70) 

L /INOUT/ TEMPS (70) , TEMPV (70) , TEMP VS (70) , QEI (70) , VX (70) , 
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S, GNDI (70) , VYl (70) , VY2 (70) ,VY3(70) ,ELSI(30), 

L VY4 (70) , VY.5 (70) , VY6 (70) , HFG , TEMPF (70) 

& /SPLIN/TIME (3) .TIMER (20) ,QSTAG(20) ,XP1 (25) , CPI (4,25) ,ND1, 

it XP2(60) , CP2 (4,60) ,ND2 ,XP3 (25) ,CP3(4,25) ,ND3,TS, 

& XP4 (25) ,CP4(4,25) ,ND4,XP5(15) ,CP5(4,15) ,ND5,XP6(15) , 

& CP6(4, 15) ,XP7 (15) ,CP7 (4,15) ,ND7,XP8 (15) ,CP8 (4,15) , 

& ND6 , ND8 , XP9 ( 15) , CP9 (4 , 15) , ND9 , XPO ( 15) , CPO (4 , 15) , NDO , 

i, ELSS(30) ,TESL(30) 

DIMENSION TA(3,96) ,QT(70) 

INTEGER TS , JK, JNT,NPP1 , FIRTS 
REWIND 8 

READ (8,*) TS , FIRTS , JK , JNT , NRES I S , IX 
READ (8,*) NPP 1 , TIMEN , RES I S , TEMPD 
IF (TS .LE. 3) THEN 
LTS - TS 
ELSE 
LTS - 3 
END IF 

READ (8 , *) (TIME (I) ,1-1, LTS) ,TVU,DEFO 
READ (8,*) ((TA(I,J) ,1-2,3) ,J-1,NP) 

IF (FIRTS .GT. 1)THEN 

READ (8,*) ( (QE (I , J) , J-l , 2) , I- 1 , NEL) 

READ (8,*) (QT (I) ,1-1, JNT) 

END IF 
RETURN 
END 


SUBROUTINE DATAOUT (TA, JK , JNT , TIMEN , 

& FIRTS, TVU.QT, TEMPD) 

£************************:M^************************************* 

C* THIS SUBROUTINE WRITES INTERMEDIATE DATA FOR RESTARTING OF * 
C* MAIN PROGRAM FROM LAST RUN * 

£**************************************************************** 


COMMON/PROC/ HRSL, ALPAR ,DIST, BI ,KREF ,ROUR, SPECR, EMIS ,HTC (140, 3) , 

£, BOLT, VOLSPR, WIDTH, QE( 140, 3) ,NRESIS,RESIS , IX 

& /ELMT/ NEL,NP ,GND(140, 3) , B (140, 3) , C (140,3) , EAREA(140) ,NPP1 , 

& IES (140,2) ,ESL (140,3) , BETA.XC (140 ,3) , YC (140, 3) , DEFO , 

it NE (140) ,TLENG, PNOD(5 , 30) ,EBC(3,3) ,CM(3,3) ,GNDB(70) 

& /INOUT/TEMPS (70) , TEMPV (70) , TEMPVS (70) , QEI (70) , VX (70) , 

it GNDI (70) , VYl (70) ,VY2(70) ,VY3(70) ,ELSI(30) , 

L VY4 (70) , VY5 (70) ,VY6(70) , HFG, TEMPF (70) 

4 /SPLIN/TIME (3) .TIMER (20) ,QSTAG(20) ,XP1(25) ,CP1 (4,25) ,ND1, 

it XP2 (60) ,CP2(4,60) ,ND2,XP3(25) ,CP3(4,25) ,ND3,TS, 

it XP4 (25) ,CP4(4, 25) ,ND4,XP5(15) ,CP5(4,15) ,ND5,XP6(15) , 

it CP6 (4, 15) ,XP7 (15) , CP7 (4, 15) ,ND7,XP8(15) ,CP8(4,15) , 

& ND6 ,ND8 ,XP9 (15) ,CP9(4,15) ,ND9,XP0(15) ,CP0(4,15) ,ND0, 

& ELSS(30) ,TESL(30) 

DIMENSION TA(3,96) ,QT (70) 

INTEGER TS.JK, JNT, NPP 1, FIRTS 
REWIND 8 

IF (TS .LE. 3) THEN 
LTS = TS 
ELSE 
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LTS - 3 
END IF 

WRITE (8 , *) TS , FIRTS , JK, JNT, NRESIS , IX 
WRITE (8 , *) NPP 1 , TIMEN , RES IS , TEMPD 
WRITE (8,*) (TIME (I) ,1-1, LTS) .TVU.DEFO 
WRITE (8,*) ((TA(I, J) ,1-2,3) , J-l.NP) 

IF (FIRTS .GT. 1) THEN 

WRITE (8,*) ((QE(I, J) ,J«1,2) ,1-1, NEL) 
WRITE (8,*) (QT (I) ,1-1, JNT) 

END IF 
RETURN 
END 


SUBROUTINE MAINOUT (TA, TSTAR , NTS ,MN, NCOLT , NROWT) 


q * * * * * * * * * * * * * * * * * * * * * * * * * * ****** * * * * * * * * * * * * * * * * * * * * * * * * * * * -J; * * 

C * THIS SUBROUTINE PRINT OUT TEMPERATURES FOR HEAT PIPE SHELL * 

C * AND WICK * 

q ****************Vc******Vc*******************Vr**Vryf ************** 

COKMON/MATG/ AMC(140,151) ,RM(140) ,THETA(140) ,FLD(140) , 

& AM(140, 140) ,BM(140) ,CIJ(3,3) ,KIJ(3,3) 

S, /TEPS/ TEMPR, TEMPI, TMEL,DETP,DELTD,TEMPG1 (140), THETAC (140), 

& THETAF (140) ,THETAO(140) ,THETT0(140) , THETAB (70) , TEMPD 

& /PROC/ HRSL , ALPAR , DIST , B I , KREF , ROUR , SPECR , EMIS , HTC ( 1 40 , 3) , 

S, BOLT, VOLSPR, WIDTH, QE(140, 3) .NRESIS , RESTS , IX 

& /ELMT/ NEL.NP , GND (140 , 3) ,B(140,3) ,C(140,3) ,EAREA(140) ,NPP1, 

£, IES (140,2) ,ESL(140,3) ,BETA,XC(140,3) ,YC(140,3) , DEFO, 

& NE (140) , TLENG, PNOD (5 , 30) ,EBC(3,3) ,CM(3,3) ,GNDB(70) 

& /INOUT/TEMPS (70) ,TEMPV(70) ,TEMPVS(70) ,.QEI(70) ,VX(70) , 

& GND I (70) , VY1 (70) , VY2 (70) , VY3 (70) , ELSI (30) , 

& VY4 (70) , VY5 (70) , VY6 (70) ,HFG,TEMPF (70) 

& / SPLIN/TIME (3) , TIMER (20) ,QSTAG(20) ,XP1(25) ,CP1(4,25) ,ND1, 

& XP2 (60) ,CP2(4,60) ,ND2,XP3(25) ,CP3(4,25) ,ND3.,TS, 

& XP4(25) ,CP4(4,25) ,ND4,XP5(15) ,CP5(4,15) ,ND5,XP6(15) , 

L CP6 (4, 15) ,XP7 (15) ,CP7(4,15) ,ND7,XP8(15) ,CP8(4,15) , 

& ND6.ND8.XP9 (15) ,CP9(4,15) ,ND9,XP0(15) ,CP0(4,15) ,ND0, 

ELSS (30) TESL (30) 

& /OUTPT/MO , IDAY , IYEAR .TITLE , CASE , NUMBER , NO , TEMPCC , HTCC , IOUT , 

& TIMEN, DELT.DELTl.DELTP 

INTEGER TS, PNOD, IT, IOUT 
DIMENSION TA(3,96) 


C PRINT BASIC DATA ON FILE JANG 

IF (IOUT .EQ. 1) THEN 

WRITE (6, 100) MO, IDAY, I YEAR 
WRITE (6, 110) TITLE, CASE, NUMBER 
WRITE (6 , 120) NO 
WRITE (6, 130) TLENG, DIST 
WRITE (6 , 140) WIDTH, ELSI (15) 

WRITE (6 , 150) TEMPI, TSTAR 
WR ITE (6 , 160) TMEL.DETP.HRSL 
WRITE (6 , 170) TEMPCC, HTCC 
WRITE (6, 180) TEMPR , EMIS 
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WRITE (6 , 190) DELT, BETA, NTS, MN 
WRITE (6 , 200) NEL.NP 
•WRITE (6, 210) 

WRITE (6, 220) 

WRITE (6, 230) 

WRITE (6, 240) 

WRITE (6, 250) 

IOUT - 2 
END IF 

C PRINT TEMPERATURE DISTRIBUTIONS 

IF (TS .LT. 3) THEN 
IT - TS 
ELSE 
IT - 3 
END IF 

WRITE (6, 300) TIME (IT) 

WRITE (6, 3 10) 

DO 400 I - 1 , NROWT 

WRITE (6,320) (TA(IT,PNOD(I, J)) , J-l , 11) 

400 CONTINUE 

WRITE (6, 330) 

DO 420 I - 1, NROWT 

WRITE (6, 320) (TA(IT,PNOD(I f J)) , J-12.NC0LT) 

420 CONTINUE 

C PRINT TEMPERATURE DISTRIBUTIONS FOR PLOT 

TMN - TIME (IT) - TIMEN 

IF (TIME (IT) .EQ. TIMEN .OR. ABS (TMN) .LT. DELT/1000.) THEN 
WRITE (7, 340) TIME (IT) 

WRITE (7, 350) (TA(IT, J) , J-l ,NP) 

IF (TIMEN .EQ. DELT1) THEN 
TIMEN - DELTP 
ELSE 

TIMEN - TIMEN + DELTP 
END IF 
END IF 

100 FORMAT (/ //20X, 37 ( ' *' ) /20X, ' * ' , 10X, ' PROGRAM - HPMAIN' , 9X, ' *' /20X, 
& '*' ,9X, ’INPUT FILE - HPDAT' ,8X, '*' /20X, 11X, 'OUTPUT - JANG', 
i, 1 IX, '*720X, '*' ,7X, ’OUTPUT FOR PLOT - DATA' ,6X, ' * 720X, ' * ' , 

& 8X, ' FOR RESTART - RESTA' , 8X, ' /20X, 

& '*' ,8X, 'DATE : ' , 13 , ’ / ' , 13, ' / ' , 15 , 7X, /20X, 37 (' *' ) /) 

110 FORMAT (/20X,37 ('-') //20X, '**' ,6X,3a9, ’ * //20X, 37 ('-’)/) 

120 FORMAT (/ 24X, '** CASE NUMBER ',16,' **'/72('-')//) 

130 FORMAT (3X, 'TOTAL LENGTH OF HEAT PIPE (TLENG) ' , 6X, E10 . 3 , 4X, ' M ' / 

& 3X, 'HEIGHT OF VAPOR SPACE (DIST) ' , 1 IX, E10 . 3 , 4X, ' M ' ) 

140 FORMAT (3X, 'ELEMENT THICKNESS (WIDTH) ', 14X, E10 . 3 , 4X ,' M ’/ 

& 3X, 'DISTANCE BETWEEN NODES (ELSI) ' , 10X, E10. 3 , 4X, ' M'/) 

150 FORMAT (3X, 'INITIAL TEMP. (TEMPI) ', 18X.F10. 3, 4X, ' K 7 

L 3X, 'TRANSIENT TEMP. OF VAPOR (TSTAR) ' , 7X , FI 0. 3 , 4X , ' K ' ) 

160 FORMAT (3X, 'PHASE CHANGE TEMP . (TMEL) ' , 14X , F10 . 3 , 4X, ' K ' / 

L 3X, ' TEMP . DIFFERENCE FROM TMEL (DETP) ’ , 6X , F10 . 3 , 4X, ' K ' / 

4 3X,' LATENT HEAT OF PHASE CHANGE (HRSL) ' ,5X,E10. 3 , 4X, ' J/KG 7) 



170 FORMAT (3X, 'REF. TEMP. FOR CONVECTION (TEMPC) ’ ,6X, FI 0.3, AX, 'K'/ 
& 3X, 'HEAT TRANSFER COEF. (HTC) ' , 14X, F10. 3 , AX,. ' W/M**2*K' ) 

180 FORMAT (3X, 'REF. TEMP. FOR RADIATION (TEMPR) 7X, F10. 3 , AX, ’K' / 

L 3X, 'EMISSIVITY' , 28X, F10. 3 , AX/) 

190 FORMAT (3X, ’TIME STEP (DELT) ’, 23X, F10 . 3 , AX, ' SECOND'/ 

& 3X,' IMPLICIT TIME SCHEME (BETA) ’, 12X, F10 . 3 , AX/ 

& 3X, ’NUMBER OF STEP FOR IMPLICIT (NTS) ’, 6X, 110 , AX/ 

& 3X, 'TOTAL NUMBER OF TIME STEP (MN) ' , 9X, I 10 , AX/) 

200 FORMAT (3X, 'NUMBER OF ELEMENT (NEL) \ 16X.I10.4X/ 

& 3X, 'NUMBER OF NODAL POINT (NP) 13X, 110, AX) 

210 FORMAT (/3X, 'INITIAL TIME' , AOX, ' SECOND ' ) 

220 FORMAT (3X, 'FINAL TIME' ,42X, ' SECOND' ) 

230 FORMAT (/3X, 'SRU') 

2A0 FORMAT (3X, 'CPU TIME' ,4AX, ' SECOND') 

250 FORMAT (3X, 'EXPENSE' ,45X, 'DOLLARS 7/72 ('-') ////) 

300 FORMAT (2X, 83 C-')/26X, '** TIME - 1 , 2X , F9 . 3 , 2X , ' SECONDS** ' /) 
310 FORMAT (5X, ' LEADING EDGE (EVAPORATOR) ' /) 

320 FORMAT ( IX , 1 2F8 . 2) 

330 FORMAT (/55X, 'TRAILING EDGE (CONDENSER) 7) 

3 AO FORMAT (F9. 3) 

350 FORMAT (10 (F8. 2, IX)) 

RETURN 

END 
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APPENDIX C 

SAMPLE INPUT DATA 

Input data axe needed for initial and boundary conditions, dimensions of the 
heat pipe, information on nodal points and elements, properties, and operating 
conditions of the program. Input file (HPDAT) has all these data. A grid generation 
program, which was written by Dr. J. G. Hartley, was used to generate data for the 
element grid. A sample for input data file HPDAT is listed in the following pages. 
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293., 371 

,700.,!., 1.045E8,. 0044,. 0127 






1. ,10. , . 

325, 100,2, 10, 1, 

1 







0 6 

138 96 0 

1 6 

2 0 






1 2 

3 








2 2 

5 








3 2 

18 








4 -3 

18 








5 3 

5 








6 3 

3 







1 

90 86 

87 .0000 

.0000 

.0250 

.0000 


.0250 

.0020 


. 5000E+01 

.0000E+00 

.OOOOE+OO 

1 -1. 0. 

2 -1. 0. 1 

0 . 

2 0. 


2 

90 87 

91 .0000 

.0000 

.0250 

.0020 


.0000 

.0020 


.5000E+01 

. 0000E+00 

O 

o 

+ 

w 

o 

o 

o 

o 

1 -1. 0. 

2 -1. 0. 1 

0 . 

2 0. 


.3 

86 82 

83 .0250 

.0000 

.0500 

.0000 


.0500 

.0020 


.5000E+01 

.0000E+00 

. 0000E+00 

♦— » 

1 

I— * 

o 

• 

2 -1. 0. 1 

0 . 

2 0. 


4 

86 83 

87 .0250 

.0000 

.0500 

.0020 


.0250 

.0020 


. 5000E+01 

.0000E+00 

O 

O 

O 

O 

M 

+ 

O 

o 

1 -1. 0. 

2 -1. 0. 1 

0 . 

2 0. 


5 

82 78 

79 .0500 

.0000 

.0750 

. .0000 


.0750 

.0020 


.5000E+01 

.0000E+00 

.0000E+00 

1 -1. 0. 

2 -1. 0. 1 

0 . 

2 0. 


6 

82 79 

83 .0500 

.0000 

.0750 

.0020 


.0500 

.0020 


.5000E+01 

.0000E+00 

.0000E+00 

1 -1. 0. 

2 -1. 0. 1 

0 . 

2 0. 


135 

91 87 

92 .0000 

.0020 

.0250 

.0020 


.0250 

.0047 


. 4000E+01 

.0000E+00 

. 0000E+00 

1 -1. 0. 

2 -1. 0. 1 

0 . 

2 0. 


136 

91 92 

94 .0000 

.0020 

.0250 

.0047 


.0000 

.0047 


. 4000E+01 

. 0000E+00 

.0000E+00 

1 -1. 0. 

2 -1. 0. 1 

0 . 

2 0. 


137 

87 83 

88 .0250 

.0020 

.0500 

.0020 


.0500 

.0047 


. 4000E+01 

.OOOOE+OO 

.0000E+00 

1 -1. 0. 

2 -1. 0. 1 

0 . 

2 0. 


138 

87 88 

92 .0250 

.0020 

.0500 

.0047 


.0250 

.0047 


. 4000E+01 

.OOOOE+OO 

.0000E+00 

1 -1. 0. 

2 -1. 0. 1 

0 . 

2 0. 


-1 

-1 -1 -1 

-1 -1 0 . 0 . 0 . 

0 . 0 . 0 . 







90 86 82 78 74 70 
66 62 58 54 50 46 
42 38 34 30 26 22 
18 14 10 6 2 3 
0 0 0 0 0 0 


293 . 


.8 
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SCALE 

FACTOR 



21 

.0 

1 . 

.05 

.87 

.1 

.48 

.15 

.33 

.2 

.25 

.25 

.19 

.3 

.15 

.35 

.13 

.4 

.1 

.45 

.08 

.5 

.07 

.55 

.06 

.6 

.055 

.65 

.05 

.7 

.05 

.75 

.049 

.8 

.049 

.85 

.049 

.9 

.049 

.95 

.049 

1 . 

-2.6 

.049 

0 . 

0 

HEAT 

FLUX 





19 

0.0 0. 

120. 

10. E3 

240. 

30. E3 

360. 50. E3 

450. 75. E3 

600. 

100. E3 

720. 

120. E3 

840. 140. E3 

1250. 230. E3 

1300. 

236. E3 

1350. 

238. E3 

1375. 239. E3 

1400. 240. E3 

1430. 

240. E3 

1450. 

240. E3 

1500. 240. E3 

1650. 240. E3 
0.0 0.0 

1800. 

240. E3 

2000. 

240. E3 



I 
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ORIGINAL PAGE IS 


APPENDIX D 

THERMAL PROPERTIES OF SODIUM 

The properties of sodium in the vapor state are expressed in terms of temper 
atures, which axe in degree Kelvin in the following property equations: 

Saturation vapor pressure[60] [N/m 2 ]: 


P = 2.29 x 10 n x ^ x 10 ^ 
Density of sodium vapor[60] [kg/m 3 ]: 

p = 6.335 x 10 8 x x 10^ 
Viscosity of sodium vapor[61] [N-S/'m 2 ]: 


H = 6.083 x 10~ 9 xT+ 1.2606 x 10“ 5 
Enthalpy of sodium vapor[60] [J/kgj: 


h f = 271,831. - 1,595.3 x T - 0.29024 x T 2 


2.4 x 10 6 


x exp 


-13.600. 

T 


h fg = 4,636.437.26 - 180.817 x T 
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Thermal properties of sodium in the solid and liquid state[62,63’ are tabulated 


as follows: 


Density of solid sodium 


T[K] 

p[kg/m 3 ] 

273. 

972.5 

293. 

968.4 

313. 

964.2 

333. 

959.9 

353. 

955.5 

371. 

951.4 


T[K] 

p [kg/m 3 ] 

283. 

970.5 

303. 

966.3 

323. 

962.1 

343. 

957.7 

363. 

953.2 


Specific heat 


T[K] 

c [J/kg-K] 
P 

273. 

1200. 

323. 

1256. 

371. 

1364. 


Conductivity 


T[K] 

K[W/m-K] 

100. 

136. 

200. 

142. 

371. 

141. 


of solid sodium 

T[K] c p [J/kg-K] 

298. 1223. 

348. 1308. 

of solid sodium 


T[K] 

K[W/m-K] 

150. 

140. 

250. 

143. 
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Density of liquid sodium 


T[K] 

[kg/m 3 ] 

T[K] 

p [kg/m 3 ] 

373. 

927. 

473. 

904. 

573. 

882. 

673. 

859. 

773. 

834. 

873. 

809. 

973. 

783. 

1073. 

757. 



Specific heat of 

liquid sodium 


T[K] 

Cp[ J/kg-K] 

T[K] 

c p [ J/kg-K] 

371. 

1385. 

373. 

1384. 

473. 

1340. 

573. 

1305. 

673. 

1279. 

773. 

1262. 

873. 

1255. 

973. 

1255. 

1073. 

1269. 

1173. 

1289. 



Conductivity of 

liquid sodium 


T[K] 

K[W/m-K] 

T[K] 

K[W/m-K] 

473. 

81.5 

573. 

75.7 

673. 

71.0 

773. 

67.2 

873. 

63.9 

973. 

61.0 

1073. 

58.3 
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APPENDIX E 


THERMAL PROPERTIES OF HASTELLOY X 
AND STAINLESS STEEL 


Thermal properties of hastelloy x[64] for heat pipe shell and type 316 stainless 
steel [64] for wick screen are tabulated as follow. 


Conductivity of hastelloy x 


T[K] 

K[W/m-K] 

373. 

11.1 

573. 

14.7 

773. 

20.6 

973. 

22.8 


Specific heat of hastelloy x 


T[K] 

c p [ J/kg-K] 

588. 

498. 

923. 

582. 

1143. 

699. 


Conductivity of stainless steel 


T[K] 

KtW/m-K] 

80.4 

8.3 

154.8 

11.5 

247.3 

14.2 

379.1 

16.5 

516.5 

18.4 

687.8 

20.8 

855.7 

22.8 

1182.1 

26.9 


T[K] 

K[W/m-K] 

107.6 

9.7 

195.2 

12.8 

299.3 

15.2 

442.5 

17.5 

611.1 

19.8 

763.7 

21.7 

980.9 

54.7 



Specific heat of stainless steel 


T[K] 

c [J/kg-K] 
P 

71.8 

286.5 

162.9 

367.0 

270.0 

433.4 

435.1 

500.7 

736.5 

561.8 

1231.1 

648.2 


T[K] 

c [ j/kg-I 
P 

114.0 

329.7 

197.9 

392.3 

337.3 

465.0 

587.6 

540.7 

967.7 

604.5 

1468.5 

690.2 
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